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THE USE OF THE MELLIN TRANSFORM IN 
FINDING THE STRESS DISTRIBUTION IN AN 
INFINITE WEDGE 
3y C. J. TRANTER (Military College of Science, Shrivenham) 
[Received 6 October 1947] 


SUMMARY 
By making use of the Mellin transform a formal solution is obtained for the stress 
distribution in an infinite wedge under fairly general conditions of surface loading. 
The results for the particular case in which each surface is subjected to unit pressure 
for a finite distance measured from the vertex of the wedge are reduced to infinite 
integrals. These can be evaluated exactly when the wedge is a semi-infinite solid 
and are in a form suitable for numerical computation for other wedge angles. 


1. ForMAL solutions for fairly general conditions of surface loading have 
recently been obtained by Tranter and Craggs (1) and Sneddon (2) in the 
case of axially symmetrical stress. If the axis of symmetry is the z-axis 
and the loading is on the curved surfaces r = constant, Tranter and 
Craggs make use of the complex form of the Fourier integral transform, 
while for loading on the plane surfaces z = constant, Sneddon employs 
the Hankel transform. A formal solution for the stress distribution in an 
infinite wedge under similar fairly general surface tractions can be obtained 
by using the Mellin transform. Such a solution would appear to cover 
cases (e.g. discontinuities) which are excluded from the solution given by 
Timoshenko (3) for a polynomial distribution of load and to give a more 
direct result than could be obtained from an extension of Shepherd’s 
work (4) for isolated forces acting at points of the faces of the wedge. 

2. Taking cylindrical polar coordinates (r,@,z) and taking the faces of 
the wedge as 0 ta, we have to find a stress function ¢ satisfying+ 
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where f,(r), f.(r) are the normal stresses and s,(r), s,(r) the shear stresses integ 
applied to the faces of the wedge, all supposed given functions of r, 
Once ¢ has been found, the stresses op, 7,9 are given by the expressions 
shown in (2) and (3), while the third stress is given by 


0, =< += —2- (4) 
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Fic. 1. The surface loading shown is for the particular case of § 3. - 


. ° o"d ad 
Assuming ¢ is such that r?+" — (n = 0, 1, 2,3), r? — (n = 1,2), and 
er" a0" 


a3, 
é | .; , 
yp +i ~ all tend to zero as r -> 00, and writing ¢ for the Mellin transform 
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sses § jntegration by parts gives: 


f r. a lg, 1 d | 
) on 4 ( n 

ions —— P-1 dr — n = 0,1, 2), 

| ; oreo” ; . de" ; , 

; | 

td on +24, dd 

2 p-ldr +1) (n= 0,2), | 

(4) | | ” or200" r= PPT) Gon , | 
: | (6) 

P -88 

| 3° Pv ldr —p(p+1)(p+2)4, | 

j or 

; | 

Cc 1d, ' 1 1 « 

ne ‘dr = p(p+1)(p+2)(p+3)¢. | 

5 , 


Multiplying (1) by r?+, integrating with respect to r from 0 to oo and 
using (6), we find 


Ss es oth . ; 
api + Lp +2)*+ P* | 92 + p?(p+ 2)7¢ = 0. (7) 


Equations (2) and (3), after multiplication by r?+* and integration with 
respect to r from 0 to 00, give 


a , d : 

p(p+1)¢ F\(p), (p+1) w =S,(p) (@ =a), (8) 
: ‘ dd : 

p(p+1)d F,(p), (p+1) = =S,(p) (@ = —a), (9) 


where 


F ( p) ( rP+ tf (r) dr. S,(p) — 


0 0 





rP+ts (r)dr (n = 1,2). (10) 


The solution of (7) is 
¢ = Asin pé+ B cos pé+C sin(p+2)0+ D cos(p+2)0, (11) 
where A, B, C, D depend on p and «a. Substitution in (8) and (9) and 
some reduction yields 
2p(p+1)G(a, p)A 
US; (p)+S,(p)}p sin(p+-2)a—{F,(p)—F,(p)}(p+-2)cos(p+2)a, (12) 


oul —2(p+1)G(a, p)C = {8,(p)+S,(p)}sin pa—{F(p)—F,(p)}cos pa, (13) 
2p(p+1)H(a, p)B 
form = {S\(p)—S,(p)}p cos( p+ 2)a+{F\(p)+F(p)}(p+2)sin(p+2)a, (14) 
—2(p+1)H(«,p)D = {S,(p)—S,(p)}cos pa+{Fi(p)+F(p)}sin pa, (15) 
‘ where ; : : 
(5) G(a, p) (p+-1)sin 2a—sin 2(p+l1)a, | (16) 


H(«, p) J 


(p+1)sin 2a+sin 2(p+ 1)a. 
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The Mellin transform ¢ of the stress function ¢ is thus completely 
known. More interest, however, lies in the stresses themselves, and if we 
denote the Mellin transform by a ‘bar’, we have 


D } 
he P a2 a | 
(r°o9) = | pede = p(p+1)¢9, | 
0 
= of | Pd) 13, _ dd T | ” 

eh | ( or * att ne ag PY | (17) 

0 } 
_  allad , dd 

27 9) = — 2 —|rP-ldr +1] = 
(r"710) | ' al: ote (Ps 70 
0 J 


The stresses can now be found from Mellin’s inversion formula (5) and 
we have 





etic ) 
] % 7 
99 Dani | p(p+] )\or-? "dp, 
1 ¢ [ad > 
c ~t us dd : 
‘nt = 5 | (p+1 a r-P-2 dp. 
c ix J 


The formal solution is now complete, ¢ and its derivatives with respect 
to @ being found from equations (11) to (16). The line integrals in (18) can 
be evaluated in terms of infinite integrals from which numerical computa- 
tion is possible. The conversion is straightforward but somewhat tedious, 
and the process is illustrated in §3 for a particular case of the loading. 


3. The particular case considered is that in which the faces of the wedge 
are each subjected to unit pressure for a distance a measured from the 
vertex, the rest of these faces being free from normal stress and the whole 
of both faces being free from shear stress. 

For this case (10) gives 


F,(p) “ F,(p) _ ( rPtl dy —= —aqpt2 (p+2), (19) 


S,(p) = S,(p) 
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Equations (12) to (15) lead to 
A C = 0, 
p(p+1)H(a, p)B —aP**sin( p+ 2)a, (20) 
(p+2)(p+1)H(a, p)D = a?+*sin pa. 


Substitution in (11) and (18) then gives: 





, l 
sin( p+ 2)a cos pé— P sin pa.cos(p +2) 








P H(«, p)' | 

] % a\? 2 9 : Pp + 4 2 : dp | 

aad) | (| sin(p--=) ee a aes x cos(p+2)0| 7, f 

, p+2 | 

4 = = (-) | sin(p+-2)a sin pé—sin pa sin(p+ 2)6 hea | 


For values of a between 0 and 47 it is easy to show that the only zero 
of H(x,p) in the strip for which the real part of p lies between —2 and 


0 is p 1. The line integrals in (21) can be replaced by integrals from 
© to 0 and from 0 to oo along the line for which the real part of p is —1, 
less 7? times the residue at p -1. Omitting details of the algebra, we 
find 
ur 7 Sin «cos 6 : 
idee. . P(€)sin{é log(a/r)} dé, 
2a 2a-+sin 2a 
U 
. | 
7 7 sin «cos 6 r sinfé log(a/r)! 
—(a9-+-0,) : + | [P()—€Q@)] HS de— | 
2a 2a-+sin 20 1+ 
0 P (22) 
cos{é log(a/r) 
| (Qe) +EP()]° ag, | | 
1+? 
0 | 
Tr | 
-T.9 ” RE -\eos{é log(a/r)} } dé, 
a ‘ | 
0 J 
where P(é), Q(€), R(E) are given by 


(€sin 2x+-sinh 2aé) P(E) = sin(a—6)cosh(a+ 6)é+ sin(a+ 6)cosh(a—6)€, 
(€sin 2a+-sinh 2a&)Q(é) cos(a—8)sinh(a+6)é+-cos(a+6)sinh(a— 4)é, 
(€sin 2v+-sinh 2a€) R(E) = sin(a—6)sinh(«+6)é—sin(«+ @)sinh(a—@)é. | 
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The results are 
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Press, 1937), p. 7. 
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~~ sinmax dx = 


-cosma dx = 


9 9 . 
2 2ar cos 8 
= — tan } 
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For the particular case (« = $7), when the wedge becomes a semi- 
infinite solid, the integrals in (22) can be evaluated exactly by making 
use of the two results (6) 


sinh 2m 
cos 2p-+-cosh 2m’ 
(4p? > 7). 
sin 2p 


cos 2p-+-cosh 2m 





4ar cos @ r2-a2cos 20 
= r4t-2a?r2cos 20-+-a4)’ 


9 9 2-—O8 
—=|2— tana) (r <a) 


7 ied lis 


‘ (24) 


r>a), 


a r2—a? 


7 





2ar cos 6 ( asin 26 ) 
> 


r4- 2a?r2cos 20-+-a4 


which can be shown to agree with those given by Love (7), who treats 
this particular case by’an entirely different method. For other values 
of « it would appear that the stresses can only be found by evaluating 
the integrals in (22) numerically. The method developed by Filon (8) for 
trigonometric integrals of this type has proved very convenient for this 


REFERENCES 
W. Craaes, Phil. Mag. (7), 36 (1945), 241. 


1 
2. I. N. SNEppDon, Proc. Cambridge Phil. Soc. 42 (1946), 260. 

3. S. TrmosHENKO, Theory of Elasticity, Ist ed. (McGraw-Hill, 1934), p. 121. 

4. W. M. SHEPHERD, Proc. Roy. Soc. A, 148 (1935), 284. I am indebted to a referee 


5. E. C. TrrcumarssH, Introduction to the Theory of Fourier Integrals (Oxford Univ. 
6. J. Epwarps, A Treatise on the Integral Calculus (Macmillan, 1922), vol. 2, p. 276. 
7 


. A. E. H. Love, Phil. Trans. Roy. Soc. A, 228 (1929), 389. 
8. L. N. G. Firion, Proc. Roy. Soc. Edin. 49 (1928-9), 38. 


. 


A 
colu 
ce lu 
iter 
nun 
are 
are 
whi 
con 


LE 


of 
th 
co 
tic 


Ww 


qi 





emi- 


king 
Aw a] rm 1 T rT rng , 
IMPROVEMENT OF AN APPROXIMATE SET OF 
r ‘ Try b rhc T 4 ae T ATC al 
‘LATENT ROOTS AND MODAL COLUMNS OF A 
ry r T a a! ‘ r T i onl 
MATRIX BY METHODS AKIN TO THOSE OF 
‘ OTA Dry 1 nWMmMIwo r 
CLASSICAL PERTURBATION THEORY 
By H. A. JAHN (University of Birmingham) 
Received 7 October 1947] 
SUMMARY 
A method is described for simultaneously improving all the latent roots and modal 
(94 columns of a given matrix, starting from a given complete set of approximate modal 
\<4) columns. It is considered that the method will be useful as a final step in any 
iteration process of determining these quantities. The method is illustrated by a 
numerical example. The modification needed when two or more of the latent roots 
are coincident, or nearly so, is very briefly indicated. The fundamental formulae 
are akin to those of classical perturbation theory, the corresponding formulae of 
| which, for the special case of a Lagrange frequency equation, are given for 
— convenience in the Appendix. 
lues 
ting 1. The fundamental equations 
0) ‘ ; ‘ “ 
me Let a complete, linearly independent set of approximate modal vectors 
this } 
D == fat) Mee cnet «=6& = I,...,%), 
of a square matrix A of order n be given. Then operating on any 2 by 
the matrix A gives a vector Az) which can be expressed as a linear 
combination of the vectors 2(°,..., 2. If we write this linear combina- 
tion as 
ere 
Ax Ag (0) 4 t 4 ald) x9), (1) 
niv. . 
“ n. 
276. where >” denotes > , then, since the 2 are approximate modes, the Ai!) 
s s=1 


8ér 
will be first approximations to the latent roots and the a!) will be small 
quantities. 


An improvement on the original mode will then be given by 


1) 


(1) (0) 1 tye (0) (2) 
= “—— Ss AW) —)M) Mla 
T s 
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so long as none of the differences A!!)—A\? are very small (see § 6), for we 
have 


—- (1) 
, % , Ors : 4 (0) 
4 woe” 


(1) (0) 
Ax\ Ax 4 


(1) 


DO 1 Fg Op SY’ Us (AMA +S aD af 
r ~-P y= re “és 4 (1) (1) — ee hat ad ; 
3 a AMY — A i 


ee (1) a, a® q@) 
= 4S" a(14. aig SY es to 
Ly VT DD Lv NV —AD 


8 


, (1) - (1) 4 (1) 
‘ \D (0) 1 5 Ars yO) Spat . Ars ast (0) 
™\ ee 4 Ni) oe AD wt 9 
a‘? 8 :. 8 
or Axi) = AVx')+ terms of second order. 


Having found this first approximation to the latent roots and the modes, 
we may now take into account the second-order terms, writing 

(1) (2),.(1) SF (2) -A1) ‘ 

ADD = NMA + DS’ aD 2, (3) 


3s “Ss 
s 


giving A!) as the second approximation to the root, and 


(2) 


— 
(2) » wha? _s rs (1) 
x BD 4 S d (4) 


2 9)°*"Ss 
dont X\?) — A‘) 


as the second approximation to the mode. Or, in general, for the pth 
approximation 

(p-1 (p—1) , ),(p— - 

AxP—Y = APP-V YS’ a PyP-D, (5) 


8 


Nr’ 


oP) giP 1) g(P 1) (6) 
r r ) (p)"s 
La NP—AP 
for the determination of AX?) and 2”) given x?-) (r,s = 1,..., 2). 


2. Determination of the coefficients. First method 
To find the coefficients A,, a,, at any stage we have the linear vector 
equation 
A,@,+ >’ a,,%, = Az,, (7) 


or, in terms of components, 


AL Dd’ Grg%p = (Az,), (6 = 1,..., 2), (8) 


8s 


where (Az,); are the components of the vector Az,. For a given value of 
r these are n linear equations for the » unknowns 


A,,a,, (8 = 1,2,...,r—1,r+1...., 2). 
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rT we Writing 
My, 4 + Cle Unt 
D X19 Xo9 . . . ene (9) 
v1, Lon ° ° . Linn 
| the solution of these equations may be written as 
' 
yy Vy + + Mp yy (AM) Bary > + Uny 
> “12 X90 ‘ : vy-1,2 (Az,)> Ur 41,2 : ene : (10) 
“y 9 Mr—1,n (Ax )n Ur +n an 
of | as 41 (Ax) Xia ny 
fie Bee « oo Rice Ee i ig. ie 
iin D.a 12 22 1,2 ( r/2 $+1,2 no) (11) 
A ln Lon ° as vs Lan (Az,), Us+in i * Xnn 
. We see that to find the improved value of the rth latent root at any 
stage we take the matrix whose columns are the modes obtained from the pre- 
| ceding stage and replace the rth column by the components of the vector 
4) | Ax, Then the improved value of the root is the value of the deter- 
minant of this modified matrix divided by the determinant of the original 
| ° nn yo (pn) : . » . 
4] modal matrix. The coefficients a!” in the expression (6), viz. 
pul 
ian (p) 
m 1 wiP—1) 1 ’ Ars y(P—-1) 
5) } - wr TZ, XP's 
“ r s 
P for the improved modes are obtained in the same manner by replacing 
the sth column of the matrix of the modes x!?~” by the components of 
Ar,, evaluating the determinant of the matrix so modified, and dividing 
again by the determinant of the original modal matrix. 
3. The equations in matrix notation 
ctor . R 
[t may be noted that the essential computational element in the method 
" described here is the calculation of the inverse, or adjoint, of the given 
approximate modal matrix. For if 
(0) (0) | 
- a ok 
(8 7 (12) 
(0) (0) 
“in ° . " a nn 
1e Ol 


is the given approximate modal matrix and we write 


A (0) (0 
Aa = 7A Lb 
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then comparison with equation (1) shows that 


Yaw... & 
(1) (a) — 
a A a 
2 - + = 4 2 
A, | na | (14) 
(1) (1) (1) 
| om we... ss & 
i.e. A, is determined by 
XA 7 
A, = G14 = —___., (15) 
: dD 
0 
where X© is the adjoint of x and D, = |a|. 


Having found A, and hence AW (r = 1...., 
improved modal matrix x” from 


2 — gMB., (16) 


where B, is derived from A, as follows: 


n) in this manner, we find the 








(1) (1) 
1 — a3) ani 
(1) (1) 7 (1) (1) 
AD —At AY) — i 
(1) (1) 
B, = mis 1 Gnz a‘) 
= (1) )() : (i). a) ] ; ; Se 
A‘ xi AQ —Ay? |, Le. (By); = Wn w (j #1), 
(1) (1) BL). = 
Ain Ayn. 1 ( ii 
(1) G) (1) (1) = 
LAP —AY AY wf | (17) 
as is seen by comparison with equation (2). Or, in general, 
i 1 
Ax?-) = 7-04 ,, (18) 
gi?) — ap of (19) 
where B,, is obtained from A, in the same manner as B, is obtained from 
A,, ie. A@-1) 
—_ je — 
(2 FJ) (B,)i; : A-))— AP =1)? (Bi = 1. 
Ww 


This is illustrated in the following example. 
4. Application to a symmetrical third-order matrix (Elementary 
Matrices (1)+, p. 310) 


By the matrix iteration method the following roots and modes of the 
matrix A shown below have been found: 


2 2 2 0-254885 —0-95670 8-2019 
A={2 & 5], ®, = | 0584225], a, —1-29429], x, = | —5-2900], 
2 5 ll 1-0 1-0 1-0 


A, = 14-4309, A, = 26152, 


A; = 0-9539. 


+ See reference (1) at the end. 
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Suppose the following approximate modal matrix 2 had been found 
for which the adjoint X, determinant D,, and Az® are as shown: 





(14) 0-25 0-96 8-2 4-0 9-16 15:748 
7) — | 0-58 -1-3 —5:3}, %O = |—5-88 —7-95 6-081 |, 
1-0 1-0 1-0 1:88 —1-21 0-2318 
3-66 —2-52 7:8 
(15) Ax 8-4 —3°42 —5l1], 
14-4 2:58 0-9 
1 the h= |x 22-0608. Then we have for the first approximation 
xo, | 14-4308 —0-0352372 — 00608682 
y AX ss 
(16) A,= D -0-:0332898 2-61530 0-00697618 
’ + 0-00248042 —0-0000614665 0-953892 
PAY? ayy ay 
=ja) A ay, 
Lay a Ap 
é i), giving AY) = 14-4308, AY) = 2-6153, AY = 0-9539, } AY = 18-0000, for the 


roots. For the improvement of the modes, with 


f = WAM — 11-8155, AV—AD = 13-4769, AY—AW = 1-6614, 








(17) 
we find 
l 0-00298228 0-00451648 
(18) | e) = 2B, = 2] —0-00281747 l —0-00419898 
(19) | +-0-000184050 —0-0000369968 1 
from . 0-254214 -~0-959558 8-20516 
0-582687 1-29807 — §-29192 
| 0-997367 1-002945 1-0003175 
. 0-254885 0-956740 8-20256 
0-584225 1-29426 — 5-29024 
lary 11-0 1-0 1-0 
f the Proceeding to the second approximation we find 
3°99598 9-1593 15-677630 
)19 xa 5°874465 7-947675 6-1405434 |, D, = 22-047237, 
900}, | 1-878485 1-211625 0-22906397 
3°67822 2-502 7°82464 
Ax! 8430895 3:38478 —5-04608 
14-430895 2-61522 +0-95392 
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Hence 
(2 3) pl? -A20898 
A a ay xo4ro 14-430898 — - 
a® 2) aQ| = a — 2-615199 — = A,, 
a®) aQ® ye) . — — 0-953903 


giving for the second approximation to the roots 
= 14-430898, AP) = 2-615199, AP) = 0-953903, ¥ A” = 18-000000. 


5. Connexion with Rayleigh’s principle and second method 
There is another method of calculating A,, a,,(s 4 r) from the vector 
equation ’ : ; “ 
1 A,2,+- >’ a2, = At, (7) 
which brings out the connexion with Rayleigh’s principle (2) for the 
determination of latent roots given approximate modes. Taking scalar 
products with the vectors 2, (J = 1,...,n) we obtain 
(Az, %;) = A,(%,-,%,)+ >’ a,,(%_,%,) (t = 1...., 0). (20) 
8 
These, for a given value of 7, are n linear equations for the » unknowns 
A,, Gy, (8 = 1, 2,...,r—1, r+1...., n). Putting for brevity 
A, = (Az,,2)), 
6, = (z,,%), 


and writing 





(21) 


we find 
On 92 © + Opp An Oypey + + Oan 


A.A, (22) 
Ont Ong 8, r-1 ; Onr+t ° ° Onn } 
811 O19 ? 7 815 =] A, 81541 . . Orn 
| 1 i 
i i a a 
8 nt One : - On,s 1 Bos 8043 . ’ Onn 


That (22) and (23) are identical with (10) and (11) may be seen as follows. 
Since the 5,, = (x,,2,) are the scalar products of two vectors we have the 
matrix relation 


On %% + + MAD Uy UM © + Mn] PM © © Xn 
V9 . i ne (24) 
ha Be «+ « @ Rs Mea «ss 8 


n2 nn nn a ln : : Xnn 


where 2,1, %,9;-++; Zp, are the components of the vector 2,. Consequently 


A = D2, (25) 
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Similarly we have the relations 


ty, «V2 Vin] [%rn ar + + Maa (AM%,)y Bay + Mn 
l 1 Ln2 Lny M4 Len ° . Lr—10n (Ax,)» Ursin ° . Xnn 
a a) eo ae. Oa. ee 2 
ae ae a ee ee a ne es 
On} 0, 2 . : On 1 flies Pax +1 si ° Ban 
where (A2,.),, (A2,)o,-++5 (Ax,),, are the components of the vector Az,. So also 
My Ty + + Uy] [%y Fy «+ Aaa (A%)y aya - - Uy 
kn1 Ung - © Lynt Lr, Lon . - Us_1 9 (Az,),, Lsiin ° . Xan 
O14, OC - - O71 6-4 A, Orsi . Sin 
. (27) 
L9On1 =One ° ° On,s—1 y On,s+1 ° ° Onn 


Consequently equations (10) and (11) may be obtained from (22) and (23) 
by dividing by D. 

The equations (21), (22), (23) have the advantage over the corresponding 
equations (9), (10), (11), that, for a symmetrical matrix A (for which the 
exact modes are orthogonal), the matrix A will be approximately diagonal 
and also, so long as none of the latent roots are very small, the terms 
A,, = (Ax,,x,) will be larger than the non-diagonal terms A,, = (Az,, x;). 
Neglecting the smaller terms we have then from (22) 


041 -O99-+-On» A, 041+ Og9-+-Op_1 > oe ere en 
ws \ A,, = (Ax,, 2,) (28) 
Oxy (Zp) 


which is Rayleigh’s approximation for the latent roots. 


For the improvement of the modes, we have, for s < r 


> 
O14 O19 07 3-1 4451 01,541 Oy - On 
O21 9 Sos-1 Are Se,41 + Ov, Son 
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where the large terms are underlined. Retaining only the largest terms 
in the expansion of this determinant, we find 


pe oe a ee 


nn t 


i (a iy+A,(— 1/7149, 844 899++35-1,9-1 864 1,8+ 1-844 8 8 


r+1r+l**Ynn: 


Since to the first approximation we have 


A = 8,,8_...8 


“nn? 
we find, on division by A, 
A,, 4,, ay 


i=. 


sa Seg Spr Ses 


The same result is obtained for s > r. Thus in all cases 





a,, = 4145 6,,—A,, 8,5 danas 1 _| A,, 8,5 | (30) 
's N 7 
Spr 855 Or, Ses A, 8, 
To this approximation we are able to improve the roots and the modes 
without the solution of equations. 
For the iteration pro upproximation we have thus 
1) 
MP) - ; (31) 
| (p-1) (p—1) 
a‘ v A rs a (32) 
“Ts ~)a\P—-1)| A(p—1) (p—1) ite 
rr Oss AY oy 
The improved mode is g en by 
(p) 
a 
a(P) — »(p—1) , rs __—s»(p—1) 6 
ee Z NP jw" (6) 
s T 8 
which becomes, since 
A(p-1) A(p-1) ] A(p-1) o(P-D 
(p) ( = — a rr q 
AP)—P) —— 7 ___ 88 | | (33) 


"at 5) (p—1) ~ S—-1) N(p—1 (p—1) s(p—1) |? 
52D — FE—D = sP-DGP=D| 4G—v 9-0 


| A(p-1) (p—1) | /| A(p—1) (p—1) | 
xP) er y(P—1) F {| AS OP | Ay? OP \ (p-1) (34) 
oe -1) S(p—1) |/ | - —1) 8 ? 
>. || A2-Y  8@-D |/|A@-Y  ge-» || 


which gives the improvement in the modes in terms of the coefficients 





§(2-) ia (a(P-, a(P—V), (35) 
A®-D — (Ag(P-D, xP-D), (36) 


derived from the modes x?—» of the preceding approximation. 
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6. Application of the second method to a third-order matrix 
Applying these formulae to the example of §4 we have with 


0-25 0-96 8-2 3°66 —2-52 7°8 
(0) 0-58 —1:3 3 Ax = | 8-4 —3-42 —5l1], 
1-0 1-0 1-0 14-4 2-58 0-9 
11-3989 0-006 0-024 20-187 —0-0336 —0-108 
(52 ) 3-61 16 0-0] Ss ° (2 ) = —()-0336 Q- 44; 52 0-042 > 
96-33 —Q-108 0-042 91-89 
giving 


A} 14-4306, AQ) = 2-6152, AY = 0-9539, > AY = 17-9997, 


l 4(0) §(0) | 
(1) . 4149 = yi 
ays 500 50 A® 50) — 0-033277, 
l A(o) 6) | 
a‘) area! Aco) stay | = + 0°0024742, 
8fY d33| AY? o% 
l A (0) §{0) ” 
aj!) a | 35236, 
a = 39 39| AW By 
l A ‘0) §{) | : 
as} 3 l= «§=—-: 00005268, 
= = 393Q| AW 3g | 
l A) 619) 
ay) = 3(0) 9(0) rnd roe aia 10838, 
5{Y d33| Agy 953’ | 
1 | AQ 3Q| be - 
a) ; 32 32 | +4 068750. 
32 599) 5 A AH 5) 
From these we find, with 
Aw —)G 11-8154, AD —AW = 13-4767, AP —AM = 1-6613, 
the following improvement of the modes (see § 3), 
wl LOB, 
0-25 0-96 8-2 1-0 +0:0029822 ++0-00451431 
0:58 -1-3 5:3 0-0028164 1-0 —0-0041383 
1-0 1-0 1-0 _(Q-00018359 —0-00003171 1-0 
0-254209 0-959514 8-20510 
= |0°582688 1-29844 —5-29200 
0-997367 1-002950 100038 
0-254880 —0-95669 8-2020 
0-584226 —1-29462 —5-2900]- 


1-0 1-0 1-0 
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It is seen, as was to be expected, that the method described in this 8. 
paragraph converges slightly less rapidly than that employed in § 4; it 7 
has, however, the distinct advantage that the solution of the set of n linear | mo 
equations is avoided. I 

. 7 . cor 
7. Application to the Lagrange frequency equation sa 

The approximate orthogonality of the approximate modes which form | spc 
the basis of the method of §5 applies only when the matrix A is sym- 
metrical. The method is, however, easily extended to a Lagrange frequency 
»>quation of the form , 
equation of the for (—AM+K)x = 0 (37) 
or (—AI+U)x = 0, (38) 
where J is the kinetic energy matrix, A the potential energy matrix, and 

U = M-1K, (39) 
*¢ . 62.6 > . . . ° A of 
if the usual extended definition of orthogonality in terms of the kinetic 
; in 
energy is made, i.e. x,, x, are defined to be orthogonal when , 
ala 
(o., M2.) = @. (40) 
For, starting from the equation * 
A,%,+ 2 4%, = Uz, (41) For 
8 
and taking scalar products with Mzx,, we have 
A,(x,; Mzx,) | } A,(Xq Mx;) (Uz,, Mx,) (M-! Kx,., Mx) (Kx,, x). — 
(42) firs 
The formulae of § 5 will then hold, assuming again that none of the latent 
roots are very small, if we replace 
Ay = (4e,,2) by 4, = (44,,4) (43) 
and a,= (2) by H, = (¢,, Bz) = (x,,2). (44) 
K(P-v KalP-), ofp-)) J 
( onsequently XP) = - Sn = \ " i) ae. (45) 
MP (Ma-), aP-») In 


is the improved value of the root, and 


. 7(p—1) (p—1) | 7(p—1) (p-—1) 
HP) = a-D4 >") 7 ie / ae er? o 


*(p—1) (p—1) 7(p—1) (p—1) |{° § 
K > MiP K is Me 


is the improvement in the mode. These may be compared with the 
corresponding formulae of classical perturbation theory given in the 
Appendix. 
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8, The case of coincident roots 

The modification needed to the nor procedure when two or 
more latent roots are coincident, or nearly so, is outlined briefly below. 
For the sake of simplicity ae the case where the first two modes 
correspond to approximately equal roots. Let y{°, y” be the approximate 
modes corresponding to these roots, whilst 2 (r Basses n) are those corre- 


sponding to the other roots. Then we will have relations of the form 


l (0 (1) ,,(0) AL) ; z'9) a 
Ay; Mi Yi +his Ye 4 ‘J he ’ (47) 
Ay py yy? | py yf? > bss ae (48) 
n 
(0) Wi ) 
Ax = De 4 HD yO4 bY yO4 SY ald 200), (49) 
r=3 s=3 


Here the coefficients pi}, u\2, wd, wo}, AY will all be of the same order 


of magnitude, but the oe te atire a) will still be small quantities. As 


classical perturbation theory (3), before proceeding farther, we first 


(in fs 
ee 
) 


Let the transformation of y\°’, y which does this (at least to the second 


diagonalize the matrix 


order of small quantities), be given by 


0 . (0) (0) Fa 
x; C114) C12Ys’ (90) 


(OQ) 


241 


’ Cc | Cos Ye’, (51) 
ind let vA! Az 


first approximation to the first two latent roots. 


be the resulting diagonal elements and consequently the 
We shall then have 


n” 
0 1)-(0) ~ (1) (0) xo 
Ax; A\ 1 lies > af =". (52) 
$=3 
ne 
Axe AVPa+ F a xO, (53) 
s=3 
Ax AD y 0) 4 = ath) 7(0), (54) 
U 3 8 


Improved modes may then be found as before, viz. 


nt 
~ ay 


Al 0 (0) ~~ 
vy x,’ Fish , (55) 


Lay NY — AY) 
3 . 


n (1) 
‘ as. 
0 AO <s (0) xe 
Xo U5 i: , oD xD" s 9 (96) 
ae ANG 8 
oo . 


(1) 
~s a 
l (0 rs (0) RW 
H Cr? 1. =’. (57) 


r23 , a r\S ) re . 
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and the process repeated if necessary. The coefficients A,, a,, can be 


T rs 


determined by any of the methods outlined in the preceding paragraphs, 


9. Conclusions 

It is considered that the method of improving the latent roots and 
modal columns of a matrix described here will be useful as a final step 
in any iteration process of determining these quantities. The special 
method of §5 and §7 becomes invalid when one or more of the latent roots 
are very small and applies further in its present form only to a symmetrical 
matrix or to one which is the product of two symmetrical matrices. It 
should be noted that the fundamental formulae of this report are akin 
to those of classical perturbation theory (2,3), as shown in the Appendix. 
The method might in fact be described as that of perturbation theory in 
reverse, for in that theory one starts from the known latent roots and 
modal columns of a given matrix and derives those of a matrix differing 
slightly from the original matrix, whilst here one starts from approximate 
modal columns and latent roots of a given matrix and deduces improved 
values of these for the same matrix. This analogy formed the basis of 
the original derivation of the formulae given here. 
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APPENDIX 
The Corresponding Formulae of Classical Perturbation Theory 
for the Case of a Lagrangian Frequency Equation 
Let M and K be the known kinetic energy and potential energy matrices of a given 


conservative system. Let x,, A, denote the known modal columns and latent roots 


of the system, so that A, = w?, where w, is the circular frequency. Then we have 


(—A,M+K)za, = 0, (58) 
or Kz, = A. Mae, (59) 
Denoting the sealar product of two modal columns by 
(<2) > Vee Seas (60) 
i 


where 2,; are the components of the modal column or vector 2,, then the orthogonality 
with respect to the kinetic energy of the modes may be expressed by 


(Mz,, 2x.) M,Brqs (61) 


where 6,, is the usual Kronecker symbol 


8, = 0 for r#s and =1 for r=s. (62) 
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The classical perturbation problem consists in finding the changes Ad,, Az, in the 
latent roots and modes of the above system due to small changes AM, AK in the 


kinetic or potential energy matrices. Putting, for the modified system, 


AX; A,+Ada,, 
x’ Let p Cre U 3s 
we have (A, + AA,)(M+AM)+K+AK}(2,+ >’ 5 2%;) 0. 


8 


Equating the small terms of first order to zero we find 


(—A,M+K) >’ ¢,,2,+(—A, AM—MAA,+AK)a, = 0. 


Taking scalar products with z,, we obtain from this vector equation, since 


(M2254) = MypSr5» K2p5Xq) = Ap May Xq) = Ap Myr Sng = Kor Bras 


o that oy ee 
A,(AM),,—M,, AA,+(AK),, = 9. 


A(AM),, . (AK) pp 
This gives AA j-Eagrel 


M,, M,, 
from which we have 
K,,+(AK),, K,,(AM),, 
AA —-—- 
A A M,, M?, 
since A, K,,/M,,- 
On the other hand, we have 
K K (AK), K,,+(AK)>, 
M M..+(AM) j. (AM),\ 
M,,\1+- i, | 
K (AK), { (AM),,) K,,++ (AK),, Ky (4M i, 
M,., M., } M,, —_ 
Thus the modified latent root is given by 
eas K,,-+ AK, 
M,, M,,+AM, 
which may be compared with equation (45) of the text. 
+ for AA, (AM),, . (AK),, 
he equivalent form j M,, “> 


for the relative change in the latent root, derived from (70), is useful. 


(63) 
(64) 


(65) 


(66) 


The fact 


that the introduction of small coupling terms (AM),,, (AK),, between the different 
normal modes of the original system has no effect on the frequency to the first order 


of small quantities may be noted. 


To find the coefficients c,, determining the changes in the modes we take scalar 


products of the vector equation (66) with x, (s #7). This gives 
Cre( A; M,, K,.) A(AM),s T (AK),; 0, 
A(AM),,+ (AK) pg 


A, Mys - Ky, 
K,,/M, 


rr? 


> 
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giving, for the modification to the mode, 





‘(| (AK),, M),, | /| K, A 
— Bey M,, Ks M,, 
which may be compared with formula (46) of the text. 
The relation 
| K, M, \A,M. M,. 
r . M,, M,,(A,—A,) (79) 


| Kuz Mag As Meg Mug 


may be noted (compare equation (33)), so that the c¢,, are of the form 


Bice 

wen hh = A 

(compare equation (2)). the 
The case of coincident latent roots of the unmodified system and the extension whe 
of the calculations to the second approximation (second-order perturbations) will at V 
be found treated in ref. (3) or in any standard text-book of quantum theory, although subs 


the formulae are given there as applied to a complex hermitian matrix. 
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IMPROVEMENT OF APPROXIMATE LATENT 
ROOTS AND VECTORS OF A SQUARE MATRIX 
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SUMMARY 


\ proof of Jahn’s method is given in matrix notation and it is pointed out that 


the latent roots are obtained with errors of the second order of small quantities 


when the original vectors are accurate to first order only. The accuracy required 
at various stages of the computation is discussed and a slight simplification for 
subsequent approximations Is indicated. 

1. Introduction 

Iv many problems of mechanical vibration and of flutter, the accurate 
determination of frequencies and modes of vibration is an important 
matter. By the use of Rayleigh’s principle accurate frequencies may be 
found from relatively inaccurate modes; but when the modes are required 
accurately considerable labour is often necessary. In (1)+ Jahn has given 
an interesting discussion of a method for the improvement of approximate 
latent roots and the associated vectors of a square matrix, which is not 
necessarily restricted to be the matrix of a dynamical system. It appears 
to the writer that there are one or two interesting avenues of thought 


arising from the paper which deserve further consideration. 


2. Statement of Jahn’s method 

The method can be very simply stated in terms of matrices, and indeed 
the matrix presentation is indicated in $3 of the original paper. In brief, 
the problem is, given a square matrix A, which for simplicity we assume 
to have distinct latent roots,{ to find as exactly as possible the modal 
matrix « and the diagonal matrix A of the latent roots such that 

Az = zi, (1) 

starting from a given approximate matrix 2p. 

Jahn’s method is expressible as follows. Evaluate 


Xo ‘As, Ay Ta, (2) 
where A, is the matrix of the diagonal elements of x)1Az, and a is the 
matrix of the remaining elements. Since 2, is approximately equal to x, 
\, will be approximately equal to A and the quantities a,, will be small. 
See reference (1) at the end 
t It does not appear to be essential to exclude complex roots, or to assume that A is 
real. 
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Next, from a determine a matrix b such that, if A,, and A,, are the rth 
and sth (diagonal) elements of A,, 
a 
b., = —*—, 3 
a Ay —;, ( 
so that the quantities b,, are also small. When b has been found, an 
improved modal matrix can be obtained from the equation 


8 


x, = x,(L+5). (4) 
For equation (3), in matrix form, is 
a = bA,—A,), (5) 
and substitution in (2) gives 
to [Avy = A,+bA,—A,b. (6) 


Hence, from (4) and (6), 
Ax, = Ax,(I+b) 

= x(A,+b6A,—A,6)(I+)) 

= Xo(1+b)A,+2,(bA,—A,b)b 

: 2, A,+terms of order b?. (7) 

If the elements of a are of the first order of smallness, so also are those 

of b, and we have thus derived matrices x,, A,, which satisfy (1) with 
errors of the second order of smallness only. 


3. Some comments 

The presentation of §2 shows that, if one of the standard methods is 
used for the determination of x51, the improved matrix x, can be obtained 
without the evaluation of determinants. The product 271Az, is found and 
(3) used directly to find 6, when x, is found directly from (4). 

An interesting point arises from the result (7). It is that A, is correct 
to the second order of small quantities, although it was derived from (2) 
in which 2» is accurate only to the first order. The parallel with Rayleigh’s 
theorem will at once spring to mind; but this result appears to be more 
general since less restriction is placed on the nature of A. However, it 
must be emphasized that 2)! must be found accurately, i.e. so that 
29 ‘%) = J with errors of second order at most, or the result is no longer 
true: see below. 

No proof that A, is correct to the second order other than that contained 
in (7) is really necessary, but it may be helpful to give an alternative. 
Let x = x +e; then (1) yields 

Ax, = —Ae+a,A+eA 


=e. , 7 ee. enol 
Xp 1Axy = A+apteA—ay 1Ae. 
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On substitution for A on the right-hand side, 


wo 1 Ax, \+ap teA—ag (ap t+e)A(xy+e)te 
A+a9 te(A—p-Ap) 
\+ 2x9 ‘ep-"(pA—Ap), 
or %, ‘Ad, \+(pA—Ap)+<a, 'e(pA—Ap), (9) 
where p = (xp+e)1e. 


Now the diagonal terms of pA—Ap are all zero; thus the only terms 
modifying A arise from the third factor in (9) which is of the order e?. 
But if an accurate reciprocal of x, is not found, it is readily shown that 
an expression p, \—Ap will appear in (9) where p, and p, though small, 
are not equal; and first-order errors will therefore appear in the approxima- 
tion to A. If this happens, then there must, by (7), be first-order errors 
in x, also, ie. the iterative process proposed will only converge provided 
appropriate accuracy is employed at each stage. 


4, The second approximation 
It is worth remarking that in a second approximation, which requires the 
inversion of x,, it is possible to make use of the known inverse of 2); this 
is a simple numerical process involving only matrix summations and 
multiplications. Thus, for example, since 
xy x,(1+5), 
where vy is accurate to second order, 
rot = (+b) 251 
(/ b+b? os+ to I (10) 
to any desired degree of accuracy. Alternatively, since x)! as already 
found initially is an approximate reciprocal of x,, the accurate reciprocal 
can be found by the successive approximation method of §4.11 of 


Elementary Matrices (2), which converges with great rapidity. 


5. A simple example 





g 
Let A ¢ ; (11) 
8 4 
Then the exact solution of equation (1) may be written 
| 10 9 | 1-5 “hs | 1-5 al 16 |: (12) 
| 8 4][1 e -% 
i.e 
| 1-5 0-75 
l l 
g (13) 
16 0 
\= | | 
0 -2 
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Suppose the given approximate matrix 2, is 
2 -] 
ty | 


so that ‘first-order’ errors of the order of 30 per cent. have been introduced, 


(14) 





Then 
_ , 1 |] < 
Xp 3 | 9 " (15) 
and on multiplying out 
49 —5 
Fe 1 
Xo Axy i 1] a 
Hence, taking the diagonal elements 
164 0 
\; | ? 1 1? (16) 
0 “3 
so that the errors in the latent roots are relatively small. Also 
0 3 
a | 1, (17) 
3 0 
and hence, with A,—A, = 183 = 5S, 
0 £ 
b ga (18) 
1 ¢ 
ob 
- 2 —l1)fl & 101] -46 
Thus 2 | ” See my 4 
1 ifle oi} | 67 61 





When constants of proportionality are removed from the columns this 
gives approximately 


19) 
1 


in good agreement with the true answer. 


1-507 0-754 
vy ° 
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NOTES ON THE SOLUTION OF ALGEBRAIC 
LINEAR SIMULTANEOUS EQUATIONS 
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SUMMARY 


In this paper four methods of solving simultaneous equations and inverting 
es are described from the viewpoint of the practical computer. The first three 
ids can be performed on ordinary desk calculating machines; the fourth uses 

Hollerith punched card equipment. The desk methods considered are the method 
mination or ‘pivotal condensation’, the method of ‘orthogonal vectors’, and 
ethod of Choleski. The elimination method is, with slight variants, the simple 

thod taught at school, used in such a way that solutions are obtained with the 
ximum possible speed and accuracy. In the method of orthogonal vectors, 
ieable only to svmmetric matrices, a new set of variables is chosen in such 
wav that the matrix is transformed to a diagonal form, from which the solution 
juations or the inverse of the matrix is immediately obtainable. The Choleski 
od, applied here again only to symmetric matrices, expresses the square matrix 
product of two triangular matrices, the reciprocation of which is a relatively 

nple operation. This method is quicker and more accurate than the others and 

n be used. with slight modifications, in the case of unsymmetric matrices. The 
lure used with punched card equipment is essentially a mechanization of the 

nation method, with a considerable advantage in speed over the corresponding 
] thod 


1. Introduction 

Metuops for the solution of algebraic linear simultaneous equations fall 
naturally into two classes. First, there is the indirect approach, in which 
the unknowns are obtained by successive approximation. Of such methods 
the procedure of Morris (1), and the relaxation method of Southwell (2), 
ire best known and the most widely used. Secondly, there is the direct 
approach, the basic characteristic of which is the fact that any desired 
accuracy can be obtained in a single application of the process, provided 
that sufficient figures are retained throughout the computation. The best 
known method of this class is that known classically as the Gaussian 
algorithm, or more recently as the method of ‘pivotal condensation’. 
Here the matrix is reduced to a triangular form by successive elimination 
of variables, and the solution is finally obtained by a process known as 
back-substitution. Alternatively, by extensions of the elimination process, 
the matrix can be reduced to a diagonal form, and the back-substitution 
avoided. 
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There are striking differences between the two types of method quoted, 
By the direct method the solution can always be obtained (provided it 
exists), and the only limit on accuracy is the number of figures to which 
the computer is prepared or compelled to work. Indirect methods, on the 
other hand, may never converge to a solution, or may converge so slowly 
as to be impracticable. 

Again, if an indirect method converges reasonably quickly, the number 
of figures retained in the computation can be increased gradually as the 
solution improves, reducing to a minimum the labour involved. By the 
direct method, on the other hand, the precision of the final solution has 
an upper bound which is completely dependent on the number of figures 
retained at the first stage of the computation.t Finally, mistakes in the 
indirect process can merely delay the convergence to an accurate solution, 
while mistakes in direct methods may completely invalidate the solution 
obtained, so that careful checks have to be applied throughout. 

There is little doubt that indirect methods at their best are quicker and 
less laborious than direct methods when performed on desk computing 
machines. They are at their best for equations in which the diagonal 
terms are large compared with off-diagonal terms, a state of affairs 
encountered often in survey problems, statistical problems, and in the 
numerical solution of certain types of ordinary and partial differential 
equations. Unfortunately there seem to be many practical problems in 
which the equations are of a type for which indirect methods are generally 
useless and direct methods must be used. 

It is with direct methods that this paper is chiefly concerned. In section 2 
we discuss some aspects of the method of pivotal condensation in the light 
of recent experience in the Mathematics Division, National Physical 
Laboratory. In section 3 we give the theory and computational details of 
two methods applicable to symmetric matrices. The first of these seemed 
to be new, but was subsequently found to be almost identical with the 
escalator process of Morris (3). So different is the theoretical approach 
that this similarity was only noticed by comparing the numerical results 
obtained by the two methods when applied to the example given by 
Morris (3). Our presentation is accordingly given here in full detail, with 
particular emphasis on the computational layout. The second method, 
due to Choleski, is very powerful but not very well known in this country. 

In §4 a numerical example is solved by each of the methods described 
in previous sections. 

In section 5 we mention, without details, a method in which punched 


+ Though there are, of course, methods for improving the accuracy of an approximate 
inverse of a matrix, or the approximate solution of a set of linear equations. 
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card equipment is adapted for the mechanization of pivotal condensation. 
This method, several times as fast as the corresponding desk computation, 
would seem to be of considerable importance in view of the widespread 
demand for rapid methods of solving equations, calculating determinants, 
and reciprocating matrices. 

The present paper is concerned with computing methods. In a com- 
panion paper (4) a mathematical analysis is given of the rounding-off 
errors inherent in some matrix processes. The main results of this 
analysis are summarized briefly in § 6. The most striking result is that 


previous estimates of error have been far too pessimistic, and the practice 


of retaining large numbers of guarding figures has been quite unnecessary. 


2. Pivotal condensation 

A. Description of method 

Though this method is very well known and has often been described, 
the relevant information is given here for the sake of completeness, as an 
introduction to § 5, and for the further inclusion of details arising from 
recent experience in the Mathematics Division. 

Briefly, the solution of simultaneous equations requires the computa- 
tion of a column matrix x, satisfying the equation 


Ax = b, (1) 


where A is a known square matrix of order n and b a known column 
matrix. Closely allied is the problem of the reciprocation of the matrix A, 
involving the computation of a square matrix X such that 

AZ. = §. (2) 
Clearly the solution of (2) is equivalent to the solution of n equations of 
the type (1), in which b denotes successively the n columns of the unit 
matrix I, and the » solutions x then give the n columns of the reciprocal 
X of A. 

For the solution of equation (1) the process of pivotal condensation or 
elimination is generally performed in such a way that the matrix A is 
transtormed into a triangular matrix. The resulting n equations contain 
n.m—1,..., 2, 1 unknowns respectively, and are solved by an obvious 
process of ‘back-substitution’. 

For the computation of the reciprocal of a matrix (equation (2)) it is 
desirable to extend the elimination process, reducing the matrix to a 
diagonal form. This extended process increases the work of elimination 
by 50 per cent., at the expense of avoiding the back-substitution. It is 
therefore not worth while for the solution of one set of simultaneous 


equations, but undoubtedly pays for the computation of a reciprocal, or 
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when the solution of equation (1) is required for constant A and several 
different b. 

The simple elimination process proceeds by the following steps: 

(i) Write down the matrix A (nm rows and n columns) side by side with 
the matrix b (nm rows and one or more columns, according to the number 
of equations to be solved). 

(ii) Caleulate the row sums of the whole matrix (A, b), and record them 
in a separate column s, say, which will be used to serve as a check on the 
computation. 

(iii) Pick out the largest of the elements a; ; of the matrix A. Suppose 
this element, the so-called ‘pivot’, is a,,. Caleulate and record the (n—1) 


quantities 
a) 


m; — =, foralls ~k. 


Any 


(iv) Multiply the kth row (the ‘pivotal’ row) of the complete matrix 
(A, b,s) by the factors m;, and add to the corresponding ith rows. That 


is, form (1) 
a‘ a;;+M; ay; 
b\ b,;+-m;b,; | for all j and alli +k. (3) 
(1) , 
8; §; T mM; OK 


The kth row of the new matrix A® is the unchanged pivotal row of A, 
and the /th column of A has zeros in all rows except the ‘th. 

(v) Check the computation by summing the elements of the rows of 
the new matrix (A®, b™). These sums should be equal, apart from 
rounding-off errors in the last figure, to the elements in the corresponding 
rows of the transformed sum s™. 

Two courses are now possible. In the simple method, whereby the 
matrix is reduced to a triangular form, the kth row and the /th column 
are henceforth left unaltered, and are not written down again. We have 
in fact effectively eliminated the unknown 2;, leaving (n—1) equations 
in (x—1) unknowns. The elimination process is then continued, picking 
out each time the largest element of A” as pivot, until only one row and 
column, that is, a single term, remains in A®-, The n pivotal rows, 
containing respectively n, n—1.,..., 2, 1 unknowns, are then assembled 
together, and the components of the solution or solutions x are obtained 
successively by an obvious but somewhat tedious and error-provoking 
process. 

Alternatively, after the first elimination process, the kth row and Ith 
column are retained in the new matrix A® and treated in subsequent 
eliminations like any other row and column. The sequence of pivotal 
elements is exactly the same as before, the pivot at any stage being the 
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largest element in the part of the matrix excluding previous pivotal rows. 
At the final stage each row has been used once, and once only, as pivotal 
row, and in the resulting matrix A~» all the elements are zero except 
these pivotal elements, one in each row and column. Simple division of 
these elements into the corresponding elements of the final matrix bi-) 
gives immediately the solution or solutions x of the equations. This 
extension of the Gauss algorithm is usually attributed to Jordan. 

A slightly different process for avoiding the back-substitution is favoured 
by some writers (5). For the solution of equation (1), the first step would 
be to write down the coefficients in the scheme 


where I, is the unit matrix of order n. Pivots are chosen from A in 
exactly the same sequence as before, and the rows of —I,, as well as of A 
are included in the elimination process, the simple one used for the reduc- 
tion of A to a triangular form. Clearly at each stage one more row of 

I, is involved, and one more row of A is omitted so that, as in the 
previous method, the number of rows to be written down is always the 
same. The number of columns on the left of the vertical line is reduced 
by one at each stage. When the last elimination has been performed, 
nothing remains above or to the left of the vertical line, and the solution 
to the equations, or the reciprocal matrix (for b = I,,) is given in the 
bottom right-hand corner. 


An example performed by each of the three methods is given in § 4. 


B. Spe al points 

A mathematical analysis of questions concerning the error of the 
reciprocal matrix is given in a companion paper (4), summarized in § 6. 
The points discussed in this section are prompted by recent practical 
experience in the Mathematics Division of the N.P.L. 

|. The choice of the largest coefficient as pivot reduces to a minimum 
the rounding-off error by ensuring in the simple elimination process that 
the factors m, are always less than unity. In the extended processes some 
of the factors may be greater than unity. For the first method this can 
occur in the factors corresponding to previous pivotal rows, and for the 
second method in the factors corresponding to the rows of —I,. This 
leature is unavoidable, and occurs in a disguised form in the ‘substituting- 
back’ process. 


Because of the horizontal line separating (A, b) from (—I,, O) we call this, colloquially, 
e ‘below-the-ling method. 


t} 
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2. Since only a finite number of figures can be retained in the computa- 
tion, even if the initial coefficients of the matrix A are known exactly. 
the solution obtained is only an approximation. This approximation x, 


should always be checked by computing the residuals (from equation (1)) 


5 = b—Ax,. 
A more accurate solution X)+- xX, can then be obtained by finding x, from 
Ax, = 6. 


Most of the elimination has already been performed in the computation 
of Xp, so little extra effort is required. 

Somewhat similarly, a more accurate inverse can be obtained in the 
form X,+X,, where X, is an approximate solution of (2), and 


X, = X,(I—AxX,). 


This second solution is worth obtaining not only when improvement is 
required, but also in those cases when accuracy of a given order is required, 
and there is doubt about the number of correct figures in the first approxi- 
mation, of which the residuals do not always give a reliable indication. 

3. Very often the coefficients and right-hand sides of the equations are 
obtained from physical considerations and experiments, and are therefore 
known to a limited degree of accuracy. Consequently any solution, for 
which the residuals are less than the ‘tolerance’ of the original system, 
is as good as any other, and it is pointless to attempt to obtain more 
accurate solutions. For example, if the coefficients are known exactly, 
whilst the right-hand sides have an error which may be as much as 5 units 
in, say, the fourth decimal place, any solution is ‘accurate’ for which the 
residuals are less than 5 units in the fourth decimal. 

4. It seems to be customary, in solving equations by elimination, to 
keep many more extra figures than are given in the original coefficients. 
In our experience this is quite unnecessary. 

To fix ideas, consider the solution by simple elimination of a set of 
equations, working throughout with, say, six decimal places. Originally 
all the coefficients are less than unity and at least one in every row lies 
between 0-1 and 1-0. The elimination process is performed until the final 
pivot is reached. This pivot is written down to six decimal places, of 
which the last one or two figures may be incorrect, due to accumulation 
of rounding-off errors. (In spite of much that has been written to the 
contrary, it is our experience that this error is small.) 

The significance of the last pivot is this, that with well-conditioned 
equations its value is not very much less, and may even be greater, than 
the first pivot. Thus it may contain as many as six significant figures, 
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resulting in solutions correct almost to six significant figures. With ill- 
conditioned equations, this last pivot may be much smaller than the 
original—possibly the first four figures after the decimal point are all 
zeros. Consequently only a few significant figures are available in the 
final pivot, and results are certainly not obtainable to more figures. 

This does not, however, mean that more accurate solutions could be 
obtained by taking originally more than six places. For if the original 
coefticients are known accurately only to six decimals, the last pivot cannot 
be obtained to better than six decimals, and if four significant figures were 
lost in the elimination, only two remain and more are meaningless. 

[f, furthermore, the original coefficients were known only to four decimals, 
and in the working to six decimals the first four figures in the final pivot 
were zero, then the original equations have effectively no solution, being 
linearly dependent. A solution obtained with the help of more figures 
would, as far as the physical problem is concerned, be meaningless and 
even misleading. 

Our conclusion is then as follows. If the original coefficients are correct 
to n places of decimals, computations should be carried out with perhaps 
one extra figure for every ten equations, these extra figures being the 
so-called ‘guarding figures’ for the accumulation of rounding-off errors. 
Then if the first » decimals of the final pivot contains m significant figures, 
no answer to better than m significant figures can be obtained from the 
given equations to the physical problem. 

These conclusions have been borne out by experiment on a set of 
eighteen ill-conditioned equations, whose coefficients were given to five 
decimals. Several guarding figures were added, and the final pivot had 
two significant figures in its first five decimals. The coefficients were then 
rounded off to four decimals, and the elimination process repeated, only 
one significant figure being left in the final pivot. The two solutions agreed 
generally to only one significant figure, bearing out the above conclusions. 


3. Two methods applicable to symmetric matrices 


One of the disadvantages of the elimination method, even in the most 
favourable case of a symmetric matrix, is that quite a substantial part 
of the total time is taken in recording the numbers used at the various 
stages. In the two methods of this section far fewer numbers are recorded, 
and the saving in time is considerable. Improved accuracy would also 
be expected, since for every recorded figure one rounding off is performed. 
Though they apply only to symmetric matrices, any matrix can be brought 
to this form by multiplication by its transpose. 

The first method, which we have called ‘the method of orthogonal 
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vectors’, requires a word of explanation. As stated above, this seemed to us 
to be new, but in fact is almost identical with the escalator process of 
Morris (3). So different is the theoretical approach, however, that this 
similarity was only noticed by comparing actual numerical results. This 
new approach is expressible very concisely in vector and matrix notation, 
and permits of a very convenient computational layout. There is also an 
important computational difference which is explained in the next section. 


A. The method of orthogonal vectors 
In the equation Ax =b (1) 


the solution x is a column matrix or vector. The present method consists 

in expressing x as a linear combination of n column matrices or vectors 

x, in the form ” 

Mi, Tes (4) 

These vectors are chosen in a special way to satisfy the ‘orthogonality’ 

condition , ’ 4 
“AZ, = &AX,=0 ( # @), (5) 


where A is necessarily symmetric. Then from (1) and (4) we have 


Y «a, Ax, = b. (6) 
= r r 
Premultiplication by x, for all s = 1...., n, then gives the n equations 
typified by 
> a, %, Ax, = =,B. (7) 


r=1 
Equations (7) constitute a set of simultaneous equations for the quantities 
x, By equation (5), all terms in (7) vanish for s 4 r, that is, the matrix 


of coefficients is diagonal. We have in fact 


x,.b 
Ny ; , (5) 
x, AX, 
and the final solution from (4) is 
n , 
x Ss x,b x.. (9) 


r 


a x. Ax, 
The orthogonal vectors are obtained as follows. Starting with the 2 


y; 41.0.0 Q} 


{ ’ 
, = {0,0,0...., af 


we calculate in turn the vectors x,, starting with x, = y,. The vector X, 


coordinate vectors 


(10) 


is oO} 2 bv 
is given by Me = Vot A. Xp on 
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in which, since x, is to be orthogonal to x,, we use equation (5) and 


obtain : : 
x, Ay, y, Ax, “ 
de ROE ae SR. (12) 
x, AX, x, AX, 
Similarly X3 = YatAg X1+Hs Xe, 


in which, making x, orthogonal to both x, and x,, we find (since x, and 


x, are already orthogonal) 


As Ys Ax, aad Ys AX, (13) 
) x, AX, x, AX, 
Finally x, Yn Ay Xi thn Xot- 
y, Ax y,, AX, 
and A “nt Ly -—F-—, ete. (14) 
x, AX, x, AX, 


These orthogonal vectors, and hence the coefficients on the left of 
equation (7), are independent of the constant term b. The extension of 
this method to the evaluation of the reciprocal of A is therefore fairly 
easy, involving merely the computation from equation (7) of n sets of 
values of « for the n columns b of the unit matrix. The remaining calcula- 
tion follows trivially from equation (4). 

In practice, we regard this method as a device whereby the original 
matrix is transformed to an almost diagonal form. It is clearly not essential 
for equation (5) to hold: given any vectors X, we can, in theory, set up 
equations (7), solve for a, and obtain the required solution from equation 
4). In practice the orthogonality condition is used because it gives the 
simplest set of equations (7). Since we work with a finite number of 
figures, equations like (5) will not be satisfied exactly, and the off-diagonal 
terms in equation (7) may differ from zero in the number of places to 
vhich we are working. We therefore accept and record these terms, and 
solve the resulting equations (7) by iteration, which in this case is very 
rapid. For large numbers of equations, particularly when they are ill- 
onditioned, this treatment adds greatly to the accuracy attainable. Here 
ve differ from the escalator process which ignores the off-diagonal terms 
ind effectively uses equation (8) for the computation of the «’s. 

In the computations it is convenient to write column matrices like x, 
nd Ax. in rows. The matrix X whose rows are the vectors x,, is then 
seen to be a lower triangular matrix, with units in the principal diagonal. 
The numbers A. u. ete., then fit into place in the upper half of the matrix, 
\, to A, occupying the places 2 to n in row 1, ps to pw, occupying places 3 


ton in row 2, ete. In this scheme the various numbers are ideally situated 
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The orthogonality condition (5) implies that the matrix whose rows are 
the column matrices Ax, is an upper triangular matrix. The effective 
vanishing of terms below the principal diagonal is one check on the 
computation. A more powerful check on Ax,, the summation check, is 


a Az A; de ee An 
a eS 
CE 

2 1 


n 





The vectors are underlined to prevent confusion with the coefficients A, py, etc. 
b> 


obtained by calculating sx,, where s is the row formed by the sum of the 
columns of A. This quantity should be effectively equal to the sum of 
the elements of Ax,. 

The computations are performed in the following order, starting at the 
point where the rth vector x, has just been formed. 

1. Calculate x, Ax, for all s <r. This is a good check on the accuracy 
of x,, and these quantities are needed in any case, being some of the 
(practically zero) coefficients of equation (7). 

2. Calculate Ax,, checking that all elements less than the rth are 
practically zero, and carrying out also the summation check. 

3. Calculate x} Ax, (another coefficient in (7)) and the quantities of the 
type A of equation (14). (Note that y,, Ax, is merely the sth element of Ax,.) 

4. Calculate x, Ax, for all s <r. These values should also be very 
small and very nearly equal to the corresponding terms x,,Ax,. A set 
of rxr coefficients in the top left-hand corner of the matrix given by 
equation (7) has now been computed. 


5. Calculate the vector x,.,.,. 


6. Proceed in this way, and in this order, until all the vectors x, and 
all the coefficients x, Ax, have been obtained. 

7. Calculate xb for all s. 

8. Solve equation (7) for a by relaxation or iteration. 

9. Calculate the final solution x from equation (4). 


10. Calculate the residuals 
§ = b— Ax, 


and obtain if required a more accurate solution by repeating the process 


from step (7) onwards with 8 in place of b. 
When the matrix is ill-conditioned, the diagonal terms x; Ax, in the 
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equations given by (7) suffer, as did the ‘pivots’ in the elimination method, 
from loss of significant figures. In experiments with a quite ill-conditioned 
set of eighteen equations, however, a more accurate solution was obtained 
than by the elimination process, working to the same number of decimals, 
with a saving of time of about one-third. This increased accuracy is due 
to some extent to the inclusion of the small but non-zero off-diagonal 
terms in the equation (7), and the saving of time to the considerable 
decrease in the number of recorded quantities. 

An example is given in § 4. 

B. Choleski’s method 

Another method, due to Choleski, for the reciprocation of a symmetric 
matrix was pointed out to us recently by John Todd. This method is 
s0 simple in theory and straightforward in practice that its neglect in this 
country is surprising. The method consists simply in expressing the 
symmetric matrix A in the form 

A = LL’, (15) 

where L is a lower triangular matrix and L’, its transpose, an upper 


triangular matrix. Then eee (16) 


and the problem is solved. 

Methods for the reciprocation of a triangular matrix are given in 
Elementary Matrices (6). The problem merely involves the back-substitu- 
tion process mentioned in §2, for equations whose right-hand sides take 
successively the columns of the unit matrix. Further, L’-! is the transpose 
of L-!, so that only one such process is necessary. 

The computation of L is also simple. Consider the case of a third-order 
matrix. We have to find the elements «, corresponding to known elements 
1, for which the following matrix equation holds: 


i YO VU Xyy%Qy gy 43, Ayn yg 
Yo1 %e Y JL OU Mg Age] = [412 Agq Mp3] (17) 
X31 gq A%gg/ \ VU QO gg 413 Ae, gg 


Vomputationally, it is better to consider the matrix product as composed 
of the scalar products of rows of the left-hand matrix. If (r,s) denotes 
the sealar product of the rth and sth rows, we have the six independent 
equations: (1.1 


) Ay | 
(2,1) = (1,2) = ay, giving successively a,,, 1, %3), 
) 


giving successively ago, ao, 


and (3,3) = 3, ~— giViNg agp. 
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Thus only the matrix L need be recorded, and its elements are obtained 
column by column. A good check is provided by writing down an extra 
row whose elements are the sums of columns so far obtained. If this row 
is denoted by s, and the corresponding row obtained from the matrix A 
denoted by &, the identity 


nth row s’ = &,, (18) 


a 
where &,, is the nth element of X, holds for all n. As soon as a column 
is completed, and the element s, is obtained, application of the identity 
(18) for 


As pointed out by A. M. Turing, it is unnecessary to compute the 


- r gives a perfect check on the computation so far performed. 


reciprocal of L if the solution is desired to only one set of simultaneous 
equations. For equations (1) and (15) imply that 


LL’x = b. (19) 
Thus writing L’x = y, we have 
Ly = b ) 
, (20) 
L’x =y/) 


and two processes of back-substitution give the solution. 

Though this method applies only to a symmetric matrix, Turing (4) has 
shown that it can be extended, without the device of normalization, to the 
unsymmetric case in a form convenient for computation. 

For symmetric matrices, this method is undoubtedly the quickest, and 
probably the most accurate, fewer figures being recorded than in any of 
the other processes described. It has also a rather important application 
in the problem of calculating the roots A and the modes x of the matrix 
equation 

(AA—B)x = 0, (21) 
where A and B are both positive definite, symmetric matrices. By any 
method the case 


(IA—C)x = 0 (22) 


is very much more easy to deal with, and equation (21) can be transformed 
into (22) by the Choleski device. We have 


(AA—B)x = 0, 
whence (AI—A-!B)x = 0. (23) 
Now write A = LL’, A-! = L’-!L-1. Equation (23) can then be written 


(AI—L’-!L-1BL’-!L’)x = 0. (24) 
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Premultiplication by L’ then gives the equation 
(AI—L“"BL’-1)L’x = 0. (25) 


Since L-'BL’-' is symmetric, equation (25) is in the required form 
(22), the variable x having been replaced by L’x. 


4, Numerical example 

The example solved in this section is the one chosen by Morris (3) as 
an illustration of his escalator method. The complete working is given 
here for each of the three methods of elimination described in §2, and 
for the vector and Choleski methods of §3. The inverse of the matrix is 
also obtained by Choleski’s method. The computational layouts should 
be self-explanatory, but the following points may be noted: 

1. The original coefficients were given to six decimals. We have worked 
throughout with eight decimals. 

2. In the elimination methods the pivots (underlined) are always on the 
principal diagonal, so that a good deal of symmetry is maintained. This 
is always the case with a positive definite, symmetric matrix. 

3. The three elimination methods give almost exactly the same solution. 
If the substituting-back process is to be avoided, it is our opinion that the 
Jordan method is preferable to the ‘below-the-line’ method. 

4. A check on any substituting-back process, particularly desirable in 
the Choleski method, is obtainable by a method analogous to the various 
summation checks mentioned previously. Thus if the solution y has been 
obtained to the equation Ly =b, 


which occurs in the Choleski method, a good check is furnished by record- 
ing the row s, whose elements are the sums of the columns of L, and 
forming the product sy. This should then be equal to the sum of the 
elements of b. All sums recorded in the computations are required for 
some check of this kind. 

5. In the elimination method the last pivot has lost three of the original 
six significant figures. Thus in accordance with earlier remarks, if the 
coefficients are known accurately to only six figures, three figures only 
are reliable in the answers. 

6. If the original coefficients are exact, answers can be obtained to any 
degree of accuracy, and the question arises as to which of the methods 
has given the best result. Conclusive evidence is not furnished by one 
example. In particular, it is difficult to be equally ‘fair’ to all methods— 
this is not necessarily achieved by retaining the same number of decimals 
throughout. However, the accurate result has been obtained, and the 
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errors in the approximate solutions are given, for what they are worth, 
in the following table, with the results given by Morris included for 
comparison. 





Error in 5th 
decimal of | Elimination | Vector | Choleski | Morris 








| 





xy 35 13 2 23 
Xg 19 6 I 12 
X3 68 25 4 47 
M4 2 13 3 25 
Xs 44 14 2 28 
xg 25 8 I 16 





Our experience is that Choleski’s method is quicker and more accurate 
than the vector method, for symmetric matrices. Both are much superior 
to pivotal condensation, which is less accurate and involves much more 
work. For unsymmetrical matrices the vector and Choleski methods 
cannot be used in their present form, but it is likely that the modified 
Choleski method given b« ‘uring (4) is still the best both for speed and 


accuracy. The usefulne the elimination process is that it can be 
‘mechanized’ and adapt: Hollerith equipment, as described in the 


next section. 


7. In one particular the above table is unfair to the vector method. 
As already stated, we are »repared to accept non-zero off-diagonal terms 
in equations (7), that is we do not insist that quantities like x) Ax, 
should be exactly zero to e number of figures retained. We do, however, 
want to know these quan ties accurately—they are subject to rounding- 
off errors, as is seen in th computation, where x, Ax, is not always equal 
to x,Ax,. This accuracy ..n be achieved by the simple process of record- 
ing the Ax, alone to one o more extra figures (depending on the magnitude 
of the vectors x,), an almc st trivial increase in the amount of work. If one 
extra figure is kept in the present example, the errors, given in the follow- 
ing table, are seen to be very little worse than those of the Choleski method. 


Error in 5th Vector 
decimal of method 


xy 3 
Xe 2 
X3 9 
X4 4 
Xs 5 
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5. Solution by Hollerith equipment 

The present Hollerith equipment appears to be somewhat unsuitable 
for the computation of scalar products, an essential part of methods such 
as that of Choleski, and for this reason the elimination methods are more 
likely than others to be readily adaptable. 

Several such methods for the theoretical use of punched card methods 
for solving simultaneous equations have been suggested by various writers 
(refs. 8 and 9). As far as we can discover, however, these methods 
have had only limited success in practice. The theoretical mechanization 
of the elimination method is fairly straightforward, but difficulties peculiar 
to digital machines have to be overcome before the method can be said 
to be practicable. These difficulties arise through the restricted capacity 
of the Hollerith multiplying punch, and the consequent need for great 
care in fixing the position of the decimal point. This trouble is particularly 
acute in the case of ill-conditioned equations. 

The methods by which these difficulties were overcome are too lengthy 
to receive detailed treatment here. A full account, together with a detailed 
description of the complete process of solution, can be obtained from the 
authors at the National Physical Laboratory. The Hollerith equipment 
used consisted of a card punch, an eight by eight decimal multiplying 
punch, a rolling total tabulator, a reproducer summary punch, and a 
collator. A special feature of the method is the fact that tabulator and 
multiplying punch are working simultaneously throughout a large part 
of the work, effecting a considerable saving in time. 

Several sets of eighteen very ill-conditioned equations have been 
solved by this method, with a time-saving factor of about eight, com- 
pared with computation on a desk machine. 


6. Note on error analysis 

An analysis of the magnitude of the errors arising from rounding off has 
recently been carried out by A. M. Turing (4), for some methods of solving 
linear equations and inverting matrices. In all those cases where a con- 
veniently simple error estimate can be given this error does not exceed 
(order of matrix) times the error which might be caused by altering the 
coefficients in the matrix to be inverted by the maximum rounding-off 
error. This in effect means that the errors only become large when the 
inverse of the matrix has large coefficients. In the case of the Jordan 
method, for example, the maximum error £ in the inverse is given by 

E| = e|M|?n3/3, 

where ¢ is the rounding-off error in each calculation, and M is the largest 
coefficient in the inverse. 
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The statistical error is very much less than this. 
The proportionality to the square of the inverse coefficients appears 


natural when we consider the case of finding the reciprocal of a scalar: 


l l . (7) 
_-— = - ~ |-)«. 
% «+e xv(x-+e) x 


Hotelling (7) has implied that the building-up error in the elimination 


methods is so large that these methods can be used only with extreme 


caution. Turing’s analysis shows that this estimate can be reached only 


in the most exceptional circumstances, and in our practical experience on 


matrices of orders up to the twentieth, some of them very ill-conditioned, 


the errors were in fact quite small. 


The work described here has been carried out as part of the Research 


programme of the National Physical Laboratory and this paper is pub- 


lished by permission of the Director of the Laboratory. 
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SUMMARY 

The equations of motion, of continuity, and of heat-transfer by bodily motion 
are derived for a frictionless perfect-gas atmosphere on a smooth rotating spherical 
Earth of radius a. The coordinates of a point P are the vertical distance (Z) from 
P to the Earth’s surface, the great circle are (S) from a fixed point O on the Earth’s 
surface to the vertical through P, and the angle (@) between the great circle through 
O and P and the meridian of O. The equations are obtained correct to terms of 
order 1/a and so can be applied to wind-systems extending to about 1,000 km. from 
O. In the case when the motion is steady (i.e. in the operator d/dt the term in @/ét 
contributes nothing) and always parallel to the Earth’s surface, dimensionless 
variables are introduced which are the ratios of the velocity, pressure, temperature, 
etc., at any point to their values at a point on the Earth’s surface situated at the 
‘outer boundary’ of the wind-system. By considering the magnitudes of the 
constants in the resulting equations, rotary motions of the atmosphere are divided 
into ‘small-scale’, with radii up to 200 km., and ‘large-scale’ with radii up to 
1,200 km. The equations are solved by approximation for some small-scale motions 
(gradient-wind, revolving fluid, ete.) in which the heat-content of unit mass of air 
is constant, and the vertical distribution of temperature is adiabatic. For small- 
scale motions of convergence or divergence, the solution breaks down near the axis 
of rotation. For large-scale motions no vortical solutions are found, but there is 
a solution in which the wind blows perpendicularly to the meridian through 0. 
It is concluded that a possible way of escape from these difficulties is to abandon the 
assumption that the motion is always parallel to the Earth’s surface. 

1. Introduction 

In dynamical meteorology, when treating of the motion of the air from 
the equations of motion, it is usual to make two main approximations, 
one with respect to the equations of motion themselves and the other 
relative to the equation of state which defines the air as a ‘fluid’. The 
first approximation consists in regarding the earth as an infinite plane 
dise and is performed as follows: A point O in latitude 4 having been 
chosen as origin, rectangular axes Ox, Oy, Oz are set up of which Ox 1s, 
for example, directed to the east, Oy to the north, and Oz along the 
vertical. If then w is the Earth’s angular velocity about its polar axis, 





k and / are 2w cos d, 2w sin d, respectively, p and p the pressure and density | 


of the (frictionless) air, the equations of motion are, neglecting terms in @”, 


P . ; . Lép 

in ee _! oP dn +lu F 2 » dw iy = ee 

dt p 0x dt p cy dt p 02 
(1.1) 
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u,v,w) being the velocity of the fluid at the point (x,y,z) at time t. It is, 
however, rarely emphasized in the text-books (Koschmieder (1), 1933; 
V. Bjerknes e¢ al. (2), 1934; Brunt (3), 1939; Haurwitz (4), 1942) that the 
range of validity of the foregoing equations is very small indeed. For 
example, the gradient-wind equation is derived by considering the case 
when w = 0 in equations (1.1) so that the wind is parallel to the plane 
(xy. But the gradient-wind relation is applied in practice to winds parallel 
to the Earth’s surface, which can be identified with the Ozxy-plane only 
ina very restricted region round O. Taking the Earth to be a sphere of 
radius 6,370 km., points in the Ozy-plane distant 200 and 375 km., 
respectively, from O are at vertical heights above the surface of 3-1 and 
\lkm. The latter point is therefore in the tropopause if O lies in tem- 
perate latitudes. Clearly therefore the use of equations (1.1) with w = 0 
for an analysis of the horizontal motion of the air in cyclonic depressions, 
vhose radii may be anything up to 1,200 km., leaves much to be desired. 

The standard text-books usually suggest that the alternative to the use 
f the rectangular coordinates described above is the employment of 
spherical polar coordinates with pole at the centre of the Earth. Yet we 
an hardly imagine a more inappropriate set of coordinates than these 
for the treatment of tropospheric motions. We are concerned with pheno- 
mena in a fluid ‘skin’ less than 18 km. in thickness on the surface of a 
sphere 6,370 km. in radius, all observational data being given in terms 
of distances measured along the surface of the sphere. There seems there- 
fore to be little advantage in selecting a coordinate system whose pole is 
it an inaccessible point 6,370 km. from the fluid whose motion we are 
investigating! What is needed is a coordinate system in which one coordi- 
nate is measured along the Earth’s surface. In Part I we develop the 
equations of motion, the equation of continuity, and the equation of heat- 
transfer by bodily motion (conduction due to molecular or turbulent 
notion neglected) taking as the coordinates of a point P, the great circle 
stance from the origin O, the vertical passing through P and the angle 
made by the plane of the great circle through O and P with the meridian 
tO, 

We turn next to the equation of state. When solving the equations of 
motion it is too often the caset that the air is treated as if it were an 
incompressible fluid (p = constant) or even to ask if it can rotate ‘like 
i solid’. Admittedly the atmosphere is a complicated fluid to deal with 
lynamically, the heat-exchanges due to the condensation or evaporation 

+ Examples are the theory of the Rayleigh-Brunt ‘revolving fluid’ (Brunt (5), 1939), 


the wave-theory of cyclones (Haurwitz (6), 1942), or the ‘solid rotation’ sought for em- 
really by Goldie (7) (1939). 
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of water giving rise, amongst other causes, to mathematical difficulties, 
Some idealization of the problem is no doubt necessary, but to replace 
air by an incompressible fluid or, worse still, by a ‘solid’, seems to be 
going too far in the direction of simplification. An idealization closer to 
the observed nature of air is to regard it as a perfect gas of uniform 
composition in which no water-vapour is present. The equation of state 
is then p = RpT, where R refers to dry air and has the value 2-868 x 108 
cm.” sec.~? deg.-} 

The equations of Part I will therefore contain as an integral element 
the definition of air as a perfect gas. In Part II we shall find one or two 
of the simpler solutions of these equations and discuss their range of 
validity. 


2. Notation 

[t has not proved possible to select a set of symbols none of which is 
already in use for some other purpose in dynamical meteorology. We 
have employed the letter a for the Earth’s radius and b,, ¢, (r = 1, 2, 3...) 
for the constants that occur in the fundamental equations. The letters 
C and A represent constants of integration, whilst A, uw, v and A,, p,, », 
(r 1,2,3,...) denote the ‘functions of integration’ which arise in the 
solution of partial differential equations of the first order. 


Part | 
The Equations of Motion, Continuity, and Heat-transfer 


3. Coordinate systems 

The fundamental equations will eventually be expressed in terms of 
a certain set of curvilinear coordinates and the simplest method of 
deriving them is probably Lagrange’s, which has also been used by 
H. Jeffreys (8), (1942) for the tidal problem. Our first step is to define 
a number of coordinate systems and the relations between them. 

For simplicity we take the Earth to be a sphere of radius a and in Fig. |, 
C is its centre and Cx,y,2, are a set of right-handed rectangular axes 
with directions ‘fixed in space’ in Newton’s sense, Cz, being directed 
towards the North Pole. O is a point fixed in the Earth’s surface in 
latitude ¢ and longitude A, NG being the Greenwich Meridian and a the 
angle which the plane NCG makes with the plane Cx,z,. Oxyz are 
rectangular right-handed axes with Ox directed to the east, Oy to the 
north, and Oz vertical, so that Oxy is the tangent plane to the sphere at 
O. Let P be a point situated above the plane Oxy; join C to P, let CP 
intersect the Earth’s surface at E, and let the plane COE intersect the 
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plane Oxy along the line OF. Draw PF parallel to Oz and PM parallel 
to OF to meet Oz at M. Then we can describe the position of P in the 
following four ways: (i) by its rectangular coordinates (2, y,,2,) referred 
to the fixed axes with origin at C; (ii) by its rectangular coordinates 
(x,y,z) referred to the axes with origin at O of which z = FP = OM; 





(iii) by its eylindrical-polar coordinates (r,6,z) of which r = OF, z = OM, 
and @ is the angle between OF and the negative direction of the Oy-axis; 
(iv) by its ‘great circle’ coordinates (S,0,Z) of which S is the are OE, 
Zis EP, and @ is now interpreted as the angle between the meridian plane 
COD and the plane of the great circle COE. 

The relation between (2,, y,,2,) and (2, y,z) is obtained thus: If ¢ = a+, 
then the direction cosines of Ox, Oy, Oz with respect to the axes Cz,, Cy,, 
(z, are, respectively, d 

(—sin &, cos x, 0), 


(—cos ¥ sin d, —sin ¥ sin ¢, cos ¢), 
(cos us cos ¢, cos d sin &, sin ¢). 


Hence the coordinates (x,,y,,2,) of P are the projections of (x,y,z+a) 
on to Cz,, Cy,, Cz,, respectively, so that 


x, (z+a)cos ¢ cos f—y sin ¢ cos—2 sin p 
y, = (z+a)cos¢dsing—ysindsing+azcosp }. (3.1) 
z, = (z+a)sin¢d+ycos¢ 

The relation between the systems (ii) and (iii) is obviously 


x= rsin§g, y = —rcos8, 2 =m &. (3.2) 
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To express (r,6,z) in terms of (S,0,Z), we refer to Fig. 2. Here 
OF = PM=r, arcOHE=S8S, and PE=Z. Since ZECO = Sia 
PC = Z-a, and MC = z-+La, we have 


r = (Z-+a)sin(S/a), 0 = @, zta = (Z+a)cos(S/a). (3.3) | 


z 


P M 














Fic. 2. 

4. Kinetic energy of unit mass 

Suppose we have a unit mass of matter in motion at P and we denote 
its kinetic energy by K, then clearly, 

2K = #+y3+2, (4.1) 

where a dot denotes differentiation with respect to ¢. Our first step is 
to express K in terms of (#,¥,2) by using (3.1) and in doing so we note 
that, since O is fixed, ¢ is a constant whilst ~ = « = w, the (constant) 
angular velocity of the Earth about its axis. If we write 


u, = tcosf/—wx sing, Uy = &sing+wrx cosy, 
v, = yoos—wy sing, Ye = ysings+wy cosy, 
w, = 2cosf—w(z+a)sin x, W, = 2sin’+a(z+a)cos yz, 


we have 
#4, = w, cos d—?, sin d—Ug, yy = w, cos d—v, sin d+.) 
2, = ycosd+ésin ¢. 


When we introduce (4.2) into (4.1) we find, after some calculation, that 


(4.2) 


K = K,—w?*K,, (4.3) 

where 2K, = #+7°+2+ ka(z+a)—a22}—l(ay—xy), (4.4) 
k = 2w cos ¢, 1 = 2wsind, (4.5) 

2K, = —|{(z+a)cos d—y sin d}?+2?]. (4.6) 


The conversion of (4.4), (4.6) into the cylindrical-polar form follows at 
once from (3.2), and we obtain 
2K, = p24 224 224 Uf {(z +a)i—érisin 0+(z+a)r0 cos 6|+17r?6, (4.7) 


2K, = —[{(z+a)cos $+ sin ¢ cos 6}?+-r*sin?6]. (4.8) 
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Here Finally we convert (4.7), (4.8) into the ‘great circle’ coordinates by using 
Sla (3.3), which, on differentiation with respect to t, give 

* = Zsin(S/a)+(Z-+-a)(S/a)cos(S/a), 6 = 6,) (4.9) 
(3.3 Z Z cos(S/a)—(Z+a)(S a)sin(S/a). J P 


Using these formulae in (4.7), (4.8) we find 


- he ‘ e jso* «6 S : S ° 11. S Ms ie »S) 
2K, = 2+(Z+a)?\— + sin? — + k|—sin 0+} sin 2— .6 cos 6} +16 sin?— 
1 a . a a} 


} a“ ¢ 
(4.10) 
‘ : | a: ae ee ee 
2K, (Z- all cos ¢ cos— + sin d sin—cos 6} + sin?—sin? 6 (4.11) 
. a a 


a \? 


which, when substituted into (4.3), give us our final form for the kinetic 


energy of unit mass. We note that the expressions (4.10), (4.11) are exact 


and are valid for all values of S up to 2za. 


5. The equations of motion 


It is known (McConnell (9), 1931) that for a perfect fluid whose motion 


Gnove | is referred to curvilinear coordinates a” (r = 1,2,3) the equations of 
motion take the form 
(4.1 LO 
} 7 0 ‘ 5 
ep is -=F—f, (r = 1,2,3), (5.1) 
Pp CA 
note 
| where (McConnell (10), 1931) 
tant 
, d/eK\ eK 
f, = ( hen p= 1,3,3 (5.2 
dt\ ex" | ox" \ )s 
/ 
p, p are the pressure and density of the fluid, and F,, f, the components of 
external force and acceleration, respectively. In the final result, we replace 
the #” by the linear components of velocity (U,V,W) referred to our 
‘great circle’ coordinates, where 
4.2) , * ’ “9 1p , , i r me 
i S U, vig? = SO V, f= Z= W. (5.3) 
hat We put K = K,—w?K, into (5.2) and use (4.10) to obtain the following 
(4.3) formulae for UfdeJ,): 
(4.4 f d (‘ K\ oK 
(4.5) dt = i es 
(4.6 d ({(Z+a\?2 _. 1(Z+a)2 V(V : S 
— i U\_— --+-U}sin 2—+- 
vs at dt\\ a ; ¥ a S\S ° a 
(Z+a)? V S Z-+-a\., . — * 
7 + 4k (1 cos 2—)cos6+k —) W sin 0+ w?—, 
(4.7) " a S a a os 


(4.8) (5.4) 
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+- ni lw (1—cos 2 2 =) +(=") Usin27| + ot, (5.5) 
2 ) a a a ali] 
ton d [eK 
aici dt aZ 
dw en ?. 8 U lV > 
= W ria(t eae sin : +k (Tsind-+5 5sin2 cos) + 


at sin? +w? net (5.6) 
S a oZ 
These are the exact expressions for the components of acceleration; 
we proceed to obtain approximate formulae by assuming that 
(i) terms in Z/a and Z/a are negligible, le 
(ii) terms in (S/a)? are negligible. (6.7) 
Expanding sin?(S/a), sin(2S/a), cos(2.S/a) in (5.4), (5.5), (5.6) and making 
use of (5.7), we find 


7 : 

f. = aaa a) —V +w? = (5.8) 

ale + +a W—U*)cose+1(U+W =} +. wre 5.9) 

: | dt rj’ c6 

f. y k(U sin 0+V cos 6)-++-w% a. (5.10) 
C4 


The terms in w* we eliminate in the usual way by remembering that 
the Earth is a spheroid and that therefore S is, strictly speaking, the are 
of the geodesic through O, whilst Z is measured along the principal normal 
to the spheroid which passes through P. If therefore the atmosphere is 
regarded as a perfect gas moving under no forces except its pressure and 
gravity, the equations (5.1) are approximately 
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These equations, and consequently their solutions, are valid only up 

to distances from the origin for which the approximation (5.7) holds good. 


‘The usual small-scale formt of these equations given in treatises on 


dynamical meteorology is the result of still another approximation in 
which terms in l/a are neglected and Z is identified with z and S with r. 


6. The equations of continuity and of heat-transfer by bodily 

motion 

The continuity of the motion is, as usual, to be secured relative to the 
moving coordinate system (S,6,Z). The quickest way of calculating the 
equation (Michal (11), 1947) is to use the metric 

dd? = dx*?+dy?+dz?, 
which by (3.2), (3.3) can be expressed as 
d>*? = (1+ Z/a)?dS?+(Z-+a)? sin®(S/a) dé?+dZ?. 
The determinant of the coefficients of the metric is 
Z+a)*. .S . 
A= (2+0Y sin?® , (6.1) 
a* a 


and the equation of continuity is then 


iI 3 , 
op , ~* a pu’) 1 » clog A 7 
a. a ae de uv —_°-- — 0, 6.2 
at 2 ge ti DY aan 
r=1 r=1 
where r g 72 r eZ 
d : (6.3) 
uw u V/S, u W 


Hence, using (6.1), we can write (6.2) as 


Op _, ApU) | l O(pV) er +e S a) = 0 





tp AP +p(— cot — +2 
ot es S 06 eZ a 3 * Bia 


and if again we make the approximation (5.7), we obtain the approximate 


equation of continuity for our ‘great circle’ coordinates as 


Op , ApU) . 1 A(pV) , a&(pW) U 2W ( 
ik ah AE AE ie tr na oe 6.4 
co °|COUCOS S 00 ° @e] p st a . aie 


Finally, the equation of heat-transfer by bodily motion is 
J dQ = Jc,dT+pd(1/p), (6.5) 


where dQ is the quantity of heat which must be added to unit mass of 
perfect gas to alter its temperature by d7' and its specific volume by 


t i.e. the cylindrical-polar form of our eqn. (1.1). 
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d(1/p). As usual, J stands for the mechanical equivalent of heat and ¢, 
for the specific heat at constant volume. The rate of addition of heat 
during the motion is 740 _ av d 1 
“a” a ray 
and if we introduce y, the ratio of the specific heats, and write 
y—1 = R/(Je,), p= R&T, 
and transform (6.6) with the aid of the equation of continuity, we obtain 
7 ‘ » 4 » 4 OT cay T 9 4 
J dQ <= ie eU ic i eV ow 2W é 
R dt y—1 dt oS 8S S 06 eZ a 


which is our final form for the equation of heat-transfer. 


(6.6) 


(6.7) 


7. Steady horizontal motion. Dimensionless variables 
By a steady motion we shall mean one in which p, p, 7, Q, U, V, W 
are functions of (S,@,Z) only, and do not involve the time explicitly, 


whilst, in a horizontal motion, W = 0 everywhere. We do not, of course, 
suggest that the restriction W = 0 holds good for all, or even for the 


most interesting, motions found in the Earth’s atmosphere (Durst and 
Sutcliffe (12), 1938; Miles (13), 1948). We seek merely to examine how 
far we can advance in the study of the motion of a perfect-gas atmosphere 
by making this a priori assumption. 

The atmosphere is a continuous fluid and therefore a disturbance in 
any part of it must have repercussions throughout the whole. Strictly 
speaking, cyclones, tornadoes, and other wind-systems do not therefore 
have ‘boundaries’ separating them off from other wind-systems. Never- 
theless it is a fact of observation that they do have ‘boundaries’ of a kind 
which can be conceived of as geometrical surfaces, extending into the 
atmosphere, beyond which the observable effects of the systems are 
negligible. Let us therefore suppose that the origin of coordinates is 
‘inside’ the system we wish to study, and that a point L due south of the 
origin lies on the boundary of the system. We choose L to lie on the 
Earth’s surface also, so that its coordinates are 

8S = &,, d= 0, Z= 0. (7.1) 
In addition, we imagine that there is a geometrical surface parallel to that 
of the Earth and at a height A above it, and we introduce the dimensionless 
variables s = S/S, 9-8. q= BA. (7.2) 
We denote by the suffix zero the values of p, 7, U, and V at the point L 
and we define the variables o, 7, v,, and vg by 


o = p/po, ¢= TZ, v. = U Up, v9 V Vo. (7.3) 


s 


It follows that s, g and oa, 7, v,, vg all have the value unity at L. 
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We further define the following thirteen constants: 


y—l1 gh U2 ) 
by : =< mm Cy —_ 7 
y RT, RT, 
b, y—l1 kl oh — J 0 
Y RT, : RT, 
b. y—1 ay, — UV 
ze y RY RT, 
y—1 Uh IS, V, 
bh. ) ot . — Soo 74 
—y a 6 es 
r9 Y 7 | 
bs y—1 ven jn l = 0 
y RT,a RT, 
LIS, Vah kS&V, 
b; : mm Ce a ‘ml 
y RT,a RT,a 
7 7 kS§ Ug 
‘ RT, a 





If we use the perfect-gas equation in the form p = Rp)T,o7, we can 
write (5.11), (5.12), (5.13) with W = 0 as 


l o(or) ov Up OV~ MY ao ae 
_ Cy Ve —? + Cg— —2 —Cq—2 — 4 09-+-C4 8 008 8 v9, (7.5) 
o cs Cs “8s ob "2 
l (or) Cv Up CVp v.U "2 
Cy V,—- +-€,-" —°+-¢, + -$ +, v,—c,8 008 00, (7.6) 
ao sod : Os "8 @& ~ £ ii j : 
lé (oT) vb Vy is ‘ 
—- — : ee (b, v, sin 06+, vg cos 0-+-b, v2 +b, v5+b; 89). 
Co Cc q Y Y 


(7.7) 
The equation of continuity (6.4) and that of heat-transfer (6.7) become 


respectively 


. log a) v Ov. - {v O(log a) 1 ov ) - 
O = tle 4 Us 1 Us|, y Je Moga) 0 8 
"\* a s' ésf' “ls 06 's af’ ait 


S, J dQ ,{ v, Ologr) , v, , 2v,) >{ v la(logr) , 1 dv 
= lt ' 8 AM I et’, |S oe =-—=<t, (7.9 
T, R dt yl ae tet tel’ "\y—is Oo 5 ae a 


In order to estimate the magnitudes of the various terms in these 
equations, we calculate the values of the constants b, and c, for the wind- 
system in a typical large depression of temperate latitudes. Goldie’s 
investigations (Goldie (14), 1939) suggest that a speed of 3 m./sec. may 
be regarded as appropriate to a point on the ‘outer boundary’ of such 
a system, whose radius may be anything up to 1,200 km. We take O to 
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be at, or near, the centre of the depression and we adopt the following 
numerical values: 


Latitude of O = 52°N., w = 7°:292« 10-5 rad./sec. 
T, = 290° A.., a = 6°37 x 108 cm. 
h = 108 em., R = 2-868 x 108 cm.? sec.-? deg.-1 


Uy = Vg = 300 cm./sec., (y—1)/y = 0-2867. 


It will be observed that we have taken the (approximate) height of the 
tropopause for h. The values of the c, are calculated for various distances 
So, whilst 6, is given for Sy = 1,200 km. only, and is therefore its maximum 
value: 

















So (km.) 1200 800 400 200 100 
10%, | 1'081 1081 1-081 1°081 1-081 
tot, =| ” ’ 
10%, | ” ” , Ps 
4 pid nacre 28 — , 
10%, =| 49°73 33°15 16°58 | 8287 | 4°144 
104¢ 5 | ” ” ’ 
1o4c, =| 7°319 3°253 0813 0°203 0-051 m 
4 > (7.10) 
10%C7 ’ ” ’ 
by = 3°383X107%h = 3-383 x 107 
b, = 9°283 x 107!2h = 9:283 x 10° 
a= » ”n » em ” 
b, = 4°871 x 104h = 4°871 x 10° 
by ’ ” ” ” 
bs = 2°239X 10713h = 2°239X 107 J 





It will be noticed that the constants c, are all small compared with 
unity, c, and c; always being the largest. For a depression of radius about 
400 km., cg, and c, are equal in magnitude to ¢,, cy, cz, whilst even at 
200 km. they are not entirely negligible. We shall therefore define a small- 
scale motion as one whose radius does not sensibly exceed 200 km.: its 
motion is governed by equations (7.5) and (7.6) with the terms in c, and 
c, omitted. A large-scale motion has a radius exceeding 200 km. and is 
controlled by equations (7.5) and (7.6) in full. 

The five ratios b,/b, (r = 1, 2,3, 4,5) are all smaller than 3 x 10-5; there- 
fore, following Jeffreys (15) (1922), we can replace (7.7) to a high degree 
of approximation by the simpler ‘hydrostatic’ equation 


sig OD a (7.11) 
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ON THE MOTION OF A PERFECT-GAS 





ATMOSPHERE 


Part II 
Some Types of Steady Horizontal Motion in the Atmosphere 


The classical investigations of Jeffreys (16, 17) (1917, 1926) and of 
G. I. Taylor (18) (1936), amongst others, have shown how to obtain 
solutions of the equations of motion of a perfect-gas atmosphere when 
the motion is unsteady and the operator d/dt reduces to 0/ét. We are here 
concerned with the other extreme case in which, in the operator d/dt, the 
term in @/ét plays no part. Our aim in solving equations (7.5), (7.6), (7.8), 
(7.9), and (7.11) is mainly illustrative to show that these equations can 
be treated by elementary mathematical methods. We shall, in the main, 
have in mind steady circular or vortex motions about a vertical axis 
through a fixed point O on the Earth’s surface. We define a symmetrical 
motion of this kind as one in which o, 7, v,, vg are functions of s and q 
alone, whilst in an asymmetrical motion, one at least of these variables 


is a function of 6 also. 


8. Simple Rotation 

The nearest analogy to a mass of air ‘spinning like a solid’ is provided 
by what we shall call ‘simple rotation’ in which v, = 0 everywhere. We 
shall prove that an adiabatic motion of this kind is possible only if (a) it 
isa small-scale motion, and (b) it is symmetrical. 

If the motion is small-scale and symmetrical, the terms in cos @ in (7.5) 
and (7.6) can be omitted, and the additional condition v, = 0 yields at 
once from (7.9) the important result dQ/dt = 0. Hence the motion is 
horizontally adiabatic. The equations (7.6), (7.8) are automatically satisfied 


whilst (7.5) reduces to 


l c(or) vB 
— = ie ae 
Oo Cs 
-_ | Cp V2 yr 
or, by (7.4), E on — +P. (8.1) 
pos S 7 


This is, of course, the gradient-wind equation in which, by hypothesis, p, 
p, and V must be functions of S and Z alone: and the heat-content of 
each unit mass is conserved in this motion. 

Suppose now that the motion is a large-scale one, the terms in cos @ in 
7.5) and (7.6) no longer being negligible. The presence of these terms 
means that we cannot now assume a@ priori that the dependent variables 
are functions of s, g alone, so that the motion is, in general, asymmetrical. 
We seek a solution in which the characteristic property of the gradient- 
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wind equation, viz. dQ/dt = 0, is still true. Equations (7.8), (7.9) then 
yield the adiabatic relation 





= silly-D (8.2) 
; i. ev 
together with Q = ——Up% nal ro, (8.3) 
’ y—1 ‘00 ali] 
whilst (7.5), (7.6) become, respectively, 
= CoP +c, Va—Cg 8 COS O v9, (8.4) 
y—lds s 
y Or Ov» . 
——_— — => —Cy Vp—. 8.5 
y—1 2 ie (8.5) 
But (8.3) and (8.5) then give 
9 Ov 
0 = (yr—c, v3)—2, 
(y7—C, V9) 20 
e Ove 2 rn . 1 °te or 
so that either > i 0 or yr = cy vg. The first possibility leads to g = ° 
C C 


which means that neither 7 nor vg can involve @ explicitly and (8.4) cannot 
then be satisfied. The second possibility implies, by (8.5), that y = —1, 
which is impossible. We therefore conclude that a large-scale simple 
rotation, analogous to the gradient wind, is not possible so long as the 
adiabatic condition for the horizontal motion is imposed. The same result 
easily follows for a small-scale asymmetrical motion: + and vg having 
been assumed to be functions of 6, as well as of s and q, the equations 
(7.8) and (7.9) can be used to show that these variables do not involve @. 

The use of the gradient-wind equation in conditions in which either the 
isobars and stream-lines are not circular or these curves are circles with 
radii greater than about 200 km., therefore implies that the heat-content 
of each unit mass of air is being altered as it moves. 


9. A quasi-adiabatic and small-scale motion of convergence or 

divergence 

The next case we want to consider is the analogue of the Rayleigh- 
Brunt (Brunt (19), 1939) rotating incompressible fluid. The motion is 
again regarded as being horizontal with a fixed vertical axis of rotation 
through O. But since v, and vg are both now different from zero, we must 
assume that along the vertical axis there is a line-sink (convergence) or 
line-source (divergence) into or out of which the necessary quantity of 
air disappears or emerges. 

We investigate whether or not there are motions in which dQ/dt = 0 
so that the equations of continuity and heat-transfer again reduce to 
a single equation together with the adiabatic relation 


c= cily-), (9.1) 
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| 
| . 
| In order to solve (7.11) we shall make the further, and independent, 
| assumption that the vertical temperature distribution corresponds to 
adiabatic equilibrium, so that we have 
T 1+ p(s, 0)—boq, (9.2) 
where (8, #) is an undetermined function of integration with the property 


that, at the point L, u(1,0) = 0. (9.3) 


We consider a small-scale symmetrical motion so that the last two 


equations become 


T 1+A,(s)—-b9q | 


: 9.4 
A,(1) f) J ( ) 


where A,(s) is an undetermined function of s alone. We can integrate 


7.8) ¢ 7.9 sive F 
7.8) and (7.9) to give osv, = Ag(q)(1—by q)!-», (9.5) 


where the function of integration has been written, for later convenience, 
4s a product in which A,(q) is arbitrary except that 


A,(0) = 1. (9.6) 
Since v, + 0, the surviving portion of (7.6) is 
ay v 
. 0 C] " 
0 C3 +¢s +C; 
cs s 
which can be integrated to cive 
A.(q) ¥ 
Vp aNd) 6, (9.7) 
Ss ° C3 


where A,(q) is another undetermined function of integration for which 


A,(0) ve. & 


(9.8) 
» c 
actin: 
Lastly, equation (7.5) becomes, in virtue of (9.1), 

‘ 9 
V OT Cy C Us v ’ 
= | 2s > 28 ie — C4 V9. 
We can transform this equation by the use of (9.4) and (9.7) and by 


observing that r 
Cs IS) 2 


2 —_ —o C4, 
C3 RT, 


to read ss / 2 a a A eae PY 


rhis equation can be integrated with respect to s and yields 


- oe A2(¢ 1c P 
/ ‘ 1 4 + 4C- a7) + = <= gt 


y—1 y/ a s* 8 Cy 


A,(q), (9.9) 
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where A,(q) is an undetermined function of q which, at L, satisfies the 








relation - . 
1 le,\? 1 ¢3) 
—A,(0) = ~{c,+c,{1+- 2 ao ott, 9 
4(9) Bat ( +52) ic,| (9.10) 
But now the equations (9.1), (9.4), and (9.5) yield 
Ao , 
v. = 2D) 1 dy q)!9-P1-+- (8) bg} 27, (9.11) 
If therefore we substitute into (9.9) we find, after some rearrangement, 
‘ 12 
<< <n Y oo C4 2 
y—1 8c 
1 6 oN" 9) ' 9 
= sales AZ(q)(1L— bp gq)? PLL -+-Aj (8) —by G} 27D 4-9 AS(q) | ++Ag(Q). 


(9.12) 


The left-hand side of this equation is a function of s alone; hence by a 


suitable choice of the hitherto undetermined functions d,(s), ,(q), A,(q), | 


and A,(q), the right-hand side must also reduce to a function of s alone. 
This is, however, impossible for real functions A,, A, because of the presence 
of the factor {1+-A,(s)—by q}-?"”-» on the right-hand side. Strictly speaking, 
there is therefore no horizontally adiabatic motion of the type we are 
seeking. But there are approximate solutions of the equation (9.12) 
obtained by using the fact that the constants c, are small compared with 
unity. We assume that the largest terms in A,(s) are of order c, and we 
neglect the squares and higher powers of these terms. We shall justify 
this approximation presently, but adopting it in the meantime, (9.12) 
becomes 
y ic. l 


——*Ay(s) — 5 8? = Ss fey Ag) +2 A5(9)} + As): 


whose right-hand side reduces to a function of s alone if 
c€, AS(q) +c, A3(q) = C?/( RT), (9.13) 
rA4(q) = —A, (9.14) 
where C? and A are (positive) constants. The function A,(s) is therefore 


y Gs 4 l ct s 


ee ee ee ee. P 9.15) 
y—1 (5) 2RT, 8s? 8c, no 
whilst, by (9.6), (9.8), and (9.10), 
C? = RT{c,+¢,(1 4-5 %) |, (9.16) 
\ ‘. 2 C5 


dims ee Ca) (9.17) 
2\RT, "4c, 
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— ys 1] o C2 | 1 (9.18) 

S the Jence, finally, “__ A,(S) 5(1—8*)| —_- — — ; 9.1% 
y—1° RT, 8? 4c, 

(9.10) fom which we see that, as we assumed earlier, A,(s) is indeed of order c,, 


since hoth C* RT. 


Our (approximate) solution for a small-scale vortical motion about a 


and c}/c, are proportional to these constants. 


vertical axis through the origin of coordinates is therefore 





(9.11) 
A,(q) l A, (s 
- p. AD) fy Als) | (9.19) 
Tb, Ss | Y l ] —byq) 
A. (¢ l Ce 
; a\f)_ *%, (9.20) 
S 2 C3 
q > | A,(s) ~bo 4s (9 21) 
(9 19 9 ) l c;\? 
I2) |} where €,A3(q) +¢2A3(9) = Cy He, r3 ); 
= Ga 
by a > 
A,(q), nd A,(s) is given by (9.18) whilst A,(q) remains arbitrary. Returning to 
alone. ir original dimensional variables (7.3) and using (7.4) we have 
sence i 9 72 , 1 Y \2 ) 
; Go? ( 0 (Vp T 519) ’ 
uking , 
7 S2— S82 12 
ve are A\(S So) = | = i"), 
» 7" 12 7 
9 19) | 2 RT, \S 
re a ‘ 7" | 
| witl 1? ( 00 (= 1 I A,(S So) | ( (9 22) 
aan ie os 7 7({? f oe 
nd we S “\h | , a 1 l- -] Z T)) 
ustify , Mo 9 ro 9/7 1 Y | 
} {C2 U3 N(Z/h)}t—4IS, 
(9.12 S . a 
T = T,+T%d,(S/S,)—TZ, J 
where [' denotes the adiabatic lapse-rate g(y—1)/(yR) and A, is an un- 
letermined function of Z/h except that d,(0) 1. To calculate the 
pressure in this vortex we write py = Rp,7) and obtain, after using (9.1) 

(¢ ) 

9.1: und (9.4), 

(9.14 p p (1 . =) . hen Y : A\(S So) | (9 23) 

0 WN a a re JL 
efore Ip i 1—] Z Ty) 
This then is the approximate solution of the equations of motion for air— 

(9.15 treated as a perfect gas—corresponding to the Rayleigh-Brunt revolving 
incompressible fluid. We may call it ‘quasi-adiabatic’ because, as we 
have seen, the adiabatic condition is not strictly satisfied. 

91 The first point to be considered is the range of validity of this solution. 
the maximum radius of the vortex is restricted to about 200 km. since, 
lor systems of greater size, the terms in c, and c, in equations (7.5) and 

0917 "2 -_ : ° 

(9.1 ‘.6) are no longer negligible. But the solution also ceases to be valid 
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in the immediate neighbourhood of the axis of rotation, for it has been 
obtained on the assumption that the square and higher powers of ),(8/S,) 
are negligible compared with unity. As an illustration, suppose that for 
cyclonic motion (Vj, > 0), we assert that our approximation begins to 


break down when 


— l y—]l 
A,(S/S,) = —— * = —0-014. 


20 vy 
? 


The corresponding value of S/S, satisfies, by (9.22), the equation 


Yo 9 9 
d a So -f poe: 1]2 
10 ay es * FS 
which, with S, = 200 km., 7, = 290° A, and the numerical values in (7.10), 
yields S = 33 km. Within this region, our approximation becomes 


increasingly inaccurate until it breaks down entirely at S = 0. This 
means that, near the axis, the motion we have found must be replaced 
by one in which either the heat-content of each unit mass is allowed to 
vary as it moves or the restriction W = 0 is removed. Difficulties associ- 
ated with the singularity at S = 0 in the solution (9.22) are therefore 
physically spurious}: the solution does not describe a possible motion of 
the air near the axis of rotation. 

The second point worthy of notice is the presence of the arbitrary 
function A,(q). By (9.5) we have 


InpSU dZ = 2p Sy Uy dZr(Z/h)\1—LZ/T,)"r-, 


so that this function controls the rate of inflow (outflow) from the line- 
sink (source) along the axis OZ. Many different kinds of vortices are 
therefore possible. For example, if A,(Z/h) > 0 for Z -Z (<h), the 
motion is a vortex with inflow (outflow) up to Z Z. which thereafter 
degenerates into the simple rotation given by 


U = 0, V = —— SB. 


Alternatively, if U) <0 and A,(Z/h) = cos(7Z/h), there is inflow up to 
Z = th and outflow beyond this height. The air is therefore sucked in to 
the OZ-axis at lower levels and ejected again at upper levels. The nature 
of A,(Z/h) cannot be determined theoretically and remains to be found 
empirically by observation of actually occurring vortices. 


Lastly we consider the march of surface-pressure and surface-tempera- 


$5, on the assumption that the solution will continue to hold up to the axis of rotation 





Compare with the conclusions reached by C. 8S. Durst and R. C. Sutcliffe (12) (1938), 
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ture from the outer edge of the vortex inwards. 


P P I + Y A (= | 
s °| y l 1 So \" 


Hence by the time we have approached to a distance of 33 km. from the 
Thus if 


is about 1,010 mb., the surface-pressure has fallen to 960 mb. Mean- 


By (9.23) the surface- 


pressure 1S 


xis of rotation the surface-pressure has fallen to p,(1—1/20). 
0 
while the surface-temperature has fallen from 7, = 290° A at the outer 


boundary to 7)(1—0-014) = 286° A at S = 33 km. The fall of pressure 
is therefore somewhat steep but is accompanied by a relatively small fall 


in temperature. 


10. A large-scale quasi-adiabatic motion 

At first sight it might seem to be a simple matter to extend the solution 
found in the last section to vortices with radii up to 1,200 km. such as 
the eyelonic depressions studied by Goldie (20) (1939). The present author 
has not succeeded in finding such solutions, and indeed the quasi-adiabatic 
motions governed by the equations (7.5) and (7.6) in full appear to be 
of a quite different character from those already found. 
table (7.1 


seven times 


Leferring to the 
0) we observe that, when S, 1,200 km., c, and c; are now some 
sreater than and c, whilst these in turn are seven times 


0 and that 


6 
ereater than c,, c,, and cz. If again we assume that dQ/dt 
the vertical temperature-distribution is adiabatic, we have 


a rity-l 11+-p(s, 0)—by gio», (10.1) 


vhilst (7.8) and (7.9) reduce to the single equation 


U;,' GSU) | y, Aor) = 
Cs oo 
i first integral of which can be written 

l Cv 
ose -(1—byq)} oe (10.2) 

U, oo 
l Woy—) OP ° 
ol (] by) ) 5,” (10.3) 


0 
vhere v is, in general, a function of s, 6, and qg. If we again assume that 
is proportional to the c,, so that its square and higher powers can be 


neo! . . . 
eglected, and if we write for brevity. 


7) ] bo: (10.4) 
e obtain , (1 ‘|S (10.5) 
¢ { er y ] 71) ral 


(10.6) 
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We now write (7.5) and (7.6) thus: 


2 
y Ou _ Ce c OV, , Cag OV’, CoH 
Pe mail gk $2080 ry + (“028+ $2 ore ‘I, 
y—l1 os | C4 C, os C4 8 oé C, 8 | 
(10.7) 
y lep Ce Co OVa . Cap Op , CaV.¥ 
ee aes — c,{v,—2scos6v, + By VO 2% Oy 3 U9 | 
y—l1s 00 | Cs C; 08 C. 8 00 Cs 8 | 
(10.8) 
and we notice that the terms enclosed in smooth brackets are multiplied 
; : Ce Cz kS, 
by constants which are equal to the square of the ratio - = — — ~ 
c Ce a 
4 5 


The fact that all the c, are still small compared with unity permits us to 
omit the terms in » when we substitute for v, and vg from (10.5), (10.6) into 
(10.7) and (10.8). We can solve the two last-mentioned equations by 
successive approximations in terms of the ratio /S,/la as follows: 


First approximation. Put v = v, and (10.7), (10.8) become 


: = Ch - Ca ~ Y_ -“ ~ Cs = (10.9) 
y—lés Wy es y—1 06 Uy ©0 
Obviously the solution of these equations is 
Y C4 Cs, ISo 
— (8,6) = Sy = vy, = =», (8, 8), (10.10) 
y—l Ny Uy RT, 


which also shows that v, is a function of s and @ alone. 


Second approximation. Put v = v,+v, and (10.7), (10.8) now read 





y Cu c, [ev Cvo C Cv 
tf “ns 
y—lés I \ és OS Cy Os , 
y lew Cz, 1fOv, , Ove Cz OV; 
dese cca: sai seein +.—_*——°3 cos 6 —-}. 
y—ls 6 U,s\cb ° 06 ¢, 06, 
The terms in v, cancel with those in pw and we are left with 
CV. Ca Cv Cv: te Cv 
—? — *scos 6 —, 2— ~scos9—. (10.11) 
Cs C4 cs C Ce C 
d o 
ey 
From these, by calculation of —— we find 
‘ cseé 
Sey Ov 
s sin 6— 14 cos 6 : =). 
cs og 
») 
whence v,(8,0) = v,(scos 0) = v(x), (10.12) 


where x = scos8@. 
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Moreover the two equations (10.11) now yield 


dv. = KSp hg 
. la dy 





10.7) kS Cr dp 
which gives us finally vo = | ~ | x ‘dy. 
: ‘ la X dy 
Third approximation. Put v = v,+v,+v, and, in the right-hand sides of 
10.8) 10.7), (10.8), substitute (i) v = v,+v2+7, in the first terms, (ii) v = v, +, 
plied in the second terms, (iii) 1 v, in the terms enclosed in smooth brackets. 
e 
hs If we also observe that 
la Cy C; ¢ i! 
"> - —— —— 7) 
us to l 0 UyN } 0 RT, 
) into we find that the terms in smooth brackets vanish identically for all functions 
gs by y, of y. It turns out that (10.7), (10.8) reduce to 
Ov. Ce CV Cy. € Ov. 
_ "3 cos @—<. ~ © ag ‘s cos d—., 
Os C4 cs fale Cs og 


and, proceeding as in the second approximation, we obtain 





10.9) 1\9 ” 
sales ( : | 2 ody, i 
“ \la] . dx 
If we repeat these operations to arrive at approximations of higher 
0.10 order we find that, at each stage, the terms in smooth brackets in (10.7) 
and (10.8) always vanish identically. Yet it is precisely these terms which 
we should expect to produce ‘vortical’ terms in the solution of the kind 
found in §9. It is in fact obvious that we shall eventually arrive at 
y—1 IS, 
L=- v1(x), 
| 7 =-_ 
i kS, \-1dv 
p= ) as. —=¢@ 
| ( la x} dx x 
Returning to our dimensional variables, the foregoing solution can be 
written as 
(S/S,)cos 8, 
0.1] x oe F ; IS, 
S ) | 
U sin (1 —" °,] “hy! > va(x) | | 
la“) dy \ > RT, a | 
, kS, \-1dv S ) 
J cos #{1——*y bs | . Io "1 \ 
la dy|  y» RT, ea} > (10.13) 
- y—1 18 - 
T=T%+4% %(x)—PZ, 
y R 
0.12 r (x) 
7 Pe Vi(x 
) (1—T'Z/T)v'9-9) 1 4- +e i 
t Po 0 i ie Ty) 
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where, since pp», 7, may be taken as the pressure and temperature at the 
point L(S». 0,0), the arbitrary function v,(x) is restricted by the condition 
v,(1) 0. 

Obviously this motion is not vortical: since U/V tan é@, the wind 
always blows roughly parallel to the great circle through O perpendicular 
to the meridian of O. The wind-speed is proportional to dv,/dy which 
itself depends on the rate of change of pressure with x. We also notice 
that, if dv,/dx is an even function of x, V changes sign as we pass through 
6 = 4m and 6 = 3x. If therefore the wind were westerly at a point 
(S,0,Z) due south of O it would be easterly at the corresponding point 
(S,7,Z) due north of O. Incidentally, we cannot jump to the conclusion 
that this solution represents the general circulation of the globe. For such 
a world-wide motion the exact expressions (5.4) to (5.6) for the acceleration 
would be required in the equations of motion instead of the approximate 
forms (5.8) to (5.10). 

Summarizing the results of §§ 9 and 10 we may say that we have obtained 
solutions of our equations under the two a priori assumptions that the 
vertical velocity-component W is everywhere zero and that the heat- 
content @ of unit mass is constant. We have seen that, under these 
conditions, the solution corresponding to small-scale vortical motion 
breaks down in the neighbourhood of the axis of rotation whilst there is 
no solution corresponding to large-scale vortical motion. But Goldie (20) 
has shown empirically that large-scale vortex motions do, as a matter of 
fact, exist in the atmosphere, his results for the surface motion of the air 
being, in our notation, 


(U2+V?)i = (constant)/S for 400 km. < S < 1,200 km. 
(U?+V?)? = (constant)S for Okm. < S < 400 km. 


Hence we conclude that, in so far as we can hope to represent atmospheric 
motions by those of a frictionless perfect gas, we must relax one or other, 
or both, of our a priori assumptions, and clearly the one which can be 
most easily dispensed with from the physical point of view is the vanishing 
of W. 
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PART I. THE GENERAL THEORY?+ 
With a Note on the Calculation of Axially-symmetrical Supersonic Flows 
by S. GOLDSTEIN 
SUMMARY 
The paper gives a review of the present state of development of the method 


of characteristics for the solution of problems of compressible flow. 


The theory applies when the flow depends on two independent variables onl 





(one of which may be the time); heat radiation, heat conduction, and viscosity are | 


neglected, and the theory in this form appears to give a good approximation in 
a large number of practical applications. The flow, if steady, is taken to b 
supersonic. 
I. Introduction 
THE method of characteristics in its present form was fully developed by 
Massau (1) in 1900. The principal mathematical ideas, and the charac- 
teristic equations in nearly the same form as that given below, will also 
be found, for example, in ref. (2). 

Until very recently (7), however, these mathematical researches re- 


° ee ‘ » ss 
mained unknown to the physicists concerned with problems of supersonic 


flow, and the theory was redeveloped ab initio. The first step was taken 
by Busemann (3), who established the well-known graphical method for 
the treatment of plane, steady, irrotational, isentropic, supersonic flow. 
His method was perfected by Temple (4), and adapted to one-dimensional 
unsteady flow by Sauer (5). However, the form and presentation of the 
theory were such that considerable difficulties were encountered in 
generalizing it either for motion with vorticity or for flow with axial 
symmetry. Progress in this direction was achieved only in recent years 
by several physicists simultaneously, the decisive step being due to 
Guderley (6), who was also the first to determine a field of flow with curved 
shock-waves. His theory has been slightly generalized, and its presenta- 
tion greatly improved, by Pfeiffer and Meyer-Koenig (7). However, it 


+ Read to the Sixth International Congress for Applied Mechanics, Paris, 1946. 
t For earlier work, cf. refs. (8) and (9). 
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view of the practical importance of the method, and because of a desire 
to stress the physical ideas to a higher degree, it was decided to attempt 
to give a new and more extensive presentation of the whole method in its 
present form. 

The use of generalized orthogonal coordinates permits a particularly 
lucid, deduction of the characteristic equations from the equations of 
motion. stressing the close connexion between the fundamental mathe- 
matical and physical properties of the characteristics. It also shows the 
essential difference between the characteristic properties of the Mach lines 
und the streamlines and leads to a simple proof, based directly on the 
physical properties of the characteristics, of the uniqueness theorem in 
its fundamental form. According to the advanced stage of generality of 
the theory this theorem, which is widely used in applications, also assumes 
1 generalized form. 

The last two sections give a similarly unified and clarified presentation 
of the cases of unsteady flow to which the method applies. Appendix II, 
by Professor 8. Goldstein, treats the role which the axis, in axially-sym- 
metrical flow, plays in the Massau method of integration of the equations 
of motion. 


II. Equations of motion for steady flow 

We shall first consider cases of steady supersonic flow. Viscosity will 
be neglected except in shock-waves, but the flow need not be vortex free. 
Heat radiation will be neglected throughout, and molecular heat conduc- 
tion also, except in shock-waves. Hence the change of state along a 
streamline will be isentropic, except where it crosses a shock-wave. 

The equation of state of the gas can be written in the form 

p = f(p.8), (1) 
where p is the pressure, p the density, and s the entropy. Hence, 
dp (ep/ép),dp+-(ép és), ds, 
where (@p/@p), = a? is the square of the speed of sound. 

It is convenient to write the equations of motion with reference to 
general orthogonal coordinates. Take any two functions a(x, y) and A(x, y) 
which are continuous and have continuous derivatives, with 

Cx OB . Ga OB 0 
Ox Ox Cy oy’ 
and denote by h(a,B)da and ha(x, B) dB 


the elements of length on the a-lines (8 = const.) and the f-lines (« = 
const.), respectively. Let «, and KB be their curvatures, taken as positive 
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when their relation to the directions of increasing « and f is as shown in The € 
Fig. 1, the direction of increasing f on the f-lines (or, briefly, the ‘normals’) 
being obtained from that of increasing « on the a-lines by a counter- 


clockwise rotation through a right angle. 
cf, re 
in th 


hadB since 











he € 
B a a 
h 
| 
Ir+h,da sine whe 
| R ='lk | | wish 
r a a | 
| | dens 
ee sen . 1. —— If 
Fic. | the 
From Fig. 1 we see that | field 
h,(«,B+dB) = h,+(eh,/e8)dB = h,(1+K, hg dp), : | a 
' 0 p 
ch,,/ep h hers 3 | 
x! Ch x Bo (3) F unk 
and similarly, Chg/ex = hy hgxg. (3) ded 
Moreover, if « denotes the angle between the a-direction and a fixed [| 10 
direction, ; 7 not 
Ky = —@e/h, @a, (4) | 
; ; : , III. 
and since the angle between the a-lines and f-lines is the same at every 
1 
point, a ‘ | 
Kg = Ce/hg op. (4) any 
a9 a9 
ro Oe O*e der 
Since a (0) F 
c as) ope X eXl 
; , | dif 
the curvatures are related to one another by 
CK CK ; 4 , 
By % ntti = 0. (6) 
h x Cx hg cp r 
Let v, and vg denote the components of the velocity in the directions " 
of increasing « and £, respectively. The condition that the entropy does 
not change along a streamline is then expressed by a2 


Cs Cs =) 


() —— = Q. (4) 
“hi Co Ph, op 
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The equation of momentum in the a-direction is 


cv ov 2 ] op 


(3) 


p h, Ox 
f. ref. 10, para. III. 39; use is also made of the equations (3)). Similarly, 


the 6-direction, 


CUpe CUe ° ] op 
v, —_ Vp = —U, Ka Vy UBKB — ° (9) 
h, Ga hg ep phgcp 
Since the radial velocity component, in axially symmetrical flow, is 
v, Sin e+ vg COS €, 
the equation ol continuity becomes 
C(pv,) , Clpvg) jp : 
wo! 4 NP 4 pv Kg t+pugk,g -/P(y sine -ugcose) = —pR, (10) 
h.éx | hg ep re 
vhere } = 0 for plane flow and j = 1 for axially symmetrical flow. (If we 


wish to allow for a distribution of sources and sinks we may include their 
lensity in the function R.) 

If the state and the velocity of the gas are given along a smooth line, L, 
the equations (1) and (7) to (10) allow us, in principle, to calculate the 
field of flow in a whole region including L. We would have to choose 
1 family of curves as our «a-lines, one of which was L. After elimination 
fp by the help of (1), the equations (7) to (10) would then determine the 
inknown normal derivatives of v,, vg, p, and s. From these we could 
deduce the values of these variables along a line L’ at a small distance 
irom L, and then continue in the same way. In general, however, this is 


+ 


not a convenient method. 


III. Mach lines 

The question which leads to a better method is the following: Are there 
iy lines on which the equations of motion say nothing about the normal 
derivatives of vg, v,, p, 8, and p? In order to show that such lines do 
exist, it is convenient to bring the equations (7) to (10) into a slightly 
lifferent form. 


Multiply (7) by (1/p)ép/és, (10) by a2/p, and add. This gives 


,f oO Ve \ ) ( ) Op sat 
q?{ a 1 OB sy eet vpn +R) +" at Bane et (7) 


h,@a 'hgop' ° ine —ph,éa' phgep- 
Multiplication of (9) by vg and subsequent subtraction from (7’) yields 
2 : a : Vp L (q2+ y2 ‘e+ a2R- UV, Cp = (v2, a?) re + V4*e], 
hy Go hy Ox ph, r hg ep 
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while multiplication of (7’) by vg and subtraction of (9), multiplied by a2 


gives 





{ev év Sa vav, Op nid 
a* vg —S — vy B + (vB +V)K UB R\ 4 3 : = —(a*—v15) M 4 
| rh, @a h, 0 p hex p he ep 
(9*) 


Multiplication of (10) by (vB —a*)/p and addition of (7*) gives 





a?" —U, UB vB + (a? v2), +a2R+-% 4 

h, a rh, @a p hy ex 

. rs 1 e(pv i o.Up ¢ 
-(vg— a2)! vgic,-+ R+ Apr.) | = (a*— v2) “P , (10* 
| p hye v | a he ep ’ 
and elimination of «g, from (7*) and (8) yields 

» s Ov, OVg “ - Pe ; v Cp | 
vs! a2 —V. Vg — + (a?#+- yy?) +a? R 4 = 

B) hy x Beat | WK p h, @a} 


9 9 ov ] Cp : 
+ (a°—vp)v, Uys u +, Ugk, 
h,éx ph, éa , 


a ov OVe , 

= (vg—a*)vgl v,—~ +vg-—+ J. (8*) 
B B “hg 8 Php op 

Evidently, the equations (7*) to (10*) do not involve any f-derivatives, | 
nor xg, if 
B Ug = +a. (11) 
Hence any region of supersonic flow is covered by two families of lines 
along which the equations of motion say nothing about the normal 
derivatives of the velocity components and the variables of state except 
for the condition contained in equation (11) itself. They are called ‘Mach 


lines’. In contrast to the ordinary lines along which, as we have seen, the 
normal derivatives are determined by the values of v,, vg, p, ¢ on the line, 
the equations of motion admit a solution with an arbitrarily chosen set 
of normal derivatives on a Mach line (subject only to (11)). On the other 
hand, when (11) is introduced into the equations (7*) to (10*), they all 
reduce to 
Vg |x a ee Vag + (r8-+ Bc 4 : am ta?R = ©. (11a) 
a a a 
which is a condition for the rate of change of the variables along the lines 
themselves. In contrast again to their properties along ordinary lines, the 
equations of motion do not admit a solution with quite arbitrary values 
of v,, Vg, p, and s along the Mach lines. 
+ It is the whole of a Mach line together with the values of the dependent variables 
along it which is called ‘characteristic’ in mathematics. In gas dynamics the term has 
been used with different meanings by different authors. 
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[II.1. In particular, a supersonic flow pattern will be compatible with 
the equations of motion which shows a discontinuity of the normal 
lerivatives of the pressure and density along a Mach line. The possibility 
{these discontinuities and their confinement to the Mach lines are the 
fundamental physical features which distinguish steady supersonic from 
steady subsonic flow. Since changes of the density gradient can be made 
visible by means of an interferometer or Schlieren-optical apparatus, the 
Mach lines can be directly observed. 

It is usual to call the Mach lines the ‘rays of propagation’ of the dis- 
ontinuities. This concept is borrowed from the theory of unsteady 
motion. In fact, it can be shown, for unsteady fluid motion without heat 
mduction and viscosity, that discontinuities of the pressure gradient are 
propagated with the local speed of sound, a, relatively to the moving 
juid.t From this it can also be deduced, for the case of steady supersonic 
low, that these discontinuities are confined to the lines given by equation 
ll). 


IV. Streamlines 


The Mach lines, however, are not the only lines along which the equa- 


tions of motion admit discontinuities of normal derivatives. If we identify 
the family of the a-lines with the streamlines, that is, if we put 


vg = Ot (12) 
in the equations (7) to (10), they reduce to 
C8 
: Q, (12a 
h Ca 
; l Aa az ¢ 
ee! Se gi A ee a (125) 
h, a ph, ex ~ he, Gx ph, ex 
by virtue of (1) and (12a)), 
. 1 op 
vn, =-—-H, (12¢) 
p hg op 
(pv,) ) ; 
(pry) | pv,kg = —IPy sine. (12d) 
h,¢ X . Fr 


These are the ‘equations of motion along the streamlines’. They show 


' For a proof ef. ref. (2), chap. VI, § 3.2. Neglecting heat conduction, heat radiation, 
hd viscosity, we are led to assume that the entropy of any individual fluid particle does 


t change with time, so that Ds/Dt 0, instead of (7). One may wish to add this equation 
the system investigated by Courant. However, it is easily seen that this does not change 
s results. 


+ It ought to be stressed that we wish to regard the equations (11) to (12d) as identities. 
dentify the family of the «-lines with the whole family of the streamlines, or the whole 
‘one of the Mach line families. The question of the regularity of the orthogonal coordinate 


‘ystems introduced in this way is investigated in Appendix I. 








202 R. E. MEYER 


that the streamlines, too, have certain characteristic properties. If we try 





to prescribe arbitrary values of v,, p, s along a streamline they will not | 
be compatible with the equations of motion unless (12a) and (125) are} 
fulfilled in addition to (12). If we have a streamline with a suitable get! 
of functions v,(«) and p(x), then xg and é@p/hg ef will be determined by 
(12c) and (12d), and hence, as is well known, there can be no jump in| 
the normal component of the pressure gradient. However, the equations 
of motion do not determine the rate of change of the velocity, v,, and the 
entropy, 8, in the direction normal to the streamlines, and discontinuities 
of these derivatives may occur.+ They correspond to the presence of 
vorticity in the flow. 





If v is the velocity vector, the equations of momentum, (8) and (9), 
can be written 
] 
grad(4\v|?)—v, curlv+-—grad p = 0. 
We introduce the enthalpy, 7, defined by 


di= Tae+™®. 


p 
Then (126) can be integrated—by virtue of (12a)—to Bernoulli’s equation 


it+4\v\? = H = const. along a streamline. (13) 


The constant, H, is- called the ‘total energy’ or ‘stagnation enthalpy’. | 


Elimination of 7 and p from these three equations gives 

vA curlv+ 7 grads—grad H = 0.t (14) 
Hence, neglecting viscosity and heat conduction, we shall find vorticity 
in the field of flow whenever the distribution of the total energy or the 








entropy is not uniform, for example, when the fluid starts from a state of | 


rest but of non-uniform temperature, or downstream of a curved shock- 
wave which gives rise to different entropy increases on different stream- 
lines. 


IV.1. The Mach lines and streamlines are, in fact, the only lines on 


which the equations of motion permit any discontinuities of normal 


derivatives. This may be seen from the equations (7*) to (10*), or by 
the following mathematical short cut which leads directly from the equa- 
tions (1) and (7) to (10) to the equations (11), (lla), and (12). If we 
consider (7) to (10) as a set of linear equations for 


ov CVea Cp Cs 
x —: —, and 5s 
hgcB hgeB hgep hg ep 

ws i } / 2] / i. | 


They are propagated with the velocity of the fluid; ef. first footnote, p. 201 
For uniform distribution of the total energy this vortex theorem was first discovered 
Dy Crocco (11). Cf. also ref. (12). 
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having expressed p by means of (1), the condition for an indeterminate 
solution will be that the rank of the matrix of this system of equations 
is less than four, that is, that all its four-row determinants vanish. This 
is easily seen to be the case if, and only if, (12) or (11) and (11a) hold. 

It is useful to note that the indeterminateness of the normal derivatives 
on the Mach lines and streamlines holds for the derivatives of any order. 
For, differentiating equations (7*) to (10*) with respect to 8 any number 
of times, we see that the f-derivatives of the highest order will always 
be multiplied by (73—a?), and those of p and v, also by Up. And it follows 
by the same kind of arguments as were used in the last paragraph that 
the Mach lines and the streamlines are the only lines on which the f- 
derivatives of any order are not fully determined by the state of motion? 


und the derivatives involving 8 to a lower order. 


Vy. Uniqueness theorem 

The boundary between a region where the state of motion is uniform 
and a region where it is not mugt be some line on which the normal 
derivatives of some order are discontinuous, and it must therefore consist 
of Mach lines and/or streamlines. Also, if two flow patterns showing no 
finite discontinuity of the velocity and the state are identical in one 
region, but not in others, this region must be bounded by Mach lines 
and/or streamlines. 

Now let us consider two flow patterns which have the same state of 
motion along a portion AB of some ordinary line, LZ. By an ordinary line 
we mean any curve in the supersonic region on which the equations of 
motion, (1) and (7) to (10), admit, and determine, one set of values of the 
normal derivatives of v,, vg, p, and s for any arbitrary (continuous) 
distribution of values of these variables themselves on the line. It follows 
irom our previous results that any curve in the supersonic region which 
does not meet any one Mach line or streamline in more than one point is 
mn ordinary line; for if it would meet, for example, a Mach line in two 
points, the equation (11@) would imply a relation between the values of 
the variables at these points. 

For both flow patterns the variables can then be expanded into the 
same Taylor series, in the neighbourhood of AB, and hence the patterns 
are identical in a small, finite region.{ Any such region, however, must 
extend to some boundary consisting of lines which are either Mach lines 
1 streamlines or sonic lines, or the axis, in a three-dimensional flow with 


By the ‘state of motion’ at a point we mean the set of the four values of the two velocity 
two independent variables of state at this point. 


> 


it the variables are analytic in a small neighbourhood of AB. 
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axial symmetry. The state of motion along AB uniquely determines the 
field of flow in the smallest ‘characteristic polygon’ which contains 4B. 

This uniqueness theorem is a more general one than is normally en- 
countered in practical cases, where the total energy and the entropy are 
known on every streamline. In that case the vorticity is determined by 





A 


ric. 2. Fic. 3a. Fic. 35. 





equation (14) and the streamlines lose their. characteristic properties. 
Except for the possible intervention of sonic lines or the axis, the unique- 
ness theorem then asserts that the state of motion on AB determines the 
flow in the ‘Mach quadrangle’ AC BD (Fig. 2). (A configuration like that 
of Fig. 3a is excluded by the definition of LZ as an ordinary line and that 
of Fig. 3b cannot represent a physically relevant flow pattern since it 
implies the intersection of different Mach lines of the same family, that is, 
two different velocities at the same point.) 


VI. Massau’s method of integration 

The important results of the previous sections, as far as the numerical 
integration of the equations of motion is concerned, may be summed up 
as follows. There are three lines through every point where the flow is 
supersonic along which the equations (11) and (12), respectively, hold. 
The state of motion on these lines cannot be chosen arbitrarily. A flow 
pattern will not be compatible with the equations of motion unless the 
equations (11a), (12a), and (125), respectively, also hold along these lines. 

The advantage of these equations lies in the fact that they are ordinary 
differential equations, their disadvantage in the fact that they refer to 
systems of coordinates which depend on the flow pattern itself. For 
numerical purposes, it is, in general, convenient to introduce Cartesian 
coordinates, x and y, in the flow plane and polar coordinates, w and 8, for 
the velocity. 

For definiteness, let the direction of increasing « on the streamlines be 
that of the velocity, and let this direction on the Mach lines always be 
chosen in such a way that it makes an acute angle with the velocity 
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direction. Equation (11) states that the angle between the streamline and 
either of the Mach lines is 
— ‘ 
p = sin-!_, (15) 
Ww 


It is called the ‘Mach angle’ (Fig. 4). If 0, «,, «_ denote the angles of the 





Fic. 4. 


velocity direction and of the tangents to the Mach lines, respectively, with 
the a-axis, the slope of the Mach lines is 


dy -_ ‘ 
tane, _ = tan(@-p).T (16) 
ax 
Since Va Ww COS LL, Vg +w sin p, Ky = —Oe/h, da, 
we have 
ov CVea 9 OlV,/V — Co +cot pone a a 
vg-—* —v, —— = vg —— ) _ w*sin*y (+ , oP ns +wdp/h, ox, 
Th, ea h, Gx * hy Ga h, x 


and equation (11a) is transformed into 


l op af pe. @ 
f : tan Ey? _ oH <u? = ll bs ad ) 
N 


+aRtanp = 0, 


p h, ¢ ¥ h,é h, ax | h, dx | 
re 
1.€, cot u = w* dé+-aRh, d« ==/Q, 
p 
Since R =1wsin 6, 
r 


and since the element of length on the Mach line is 


h,da« = dr/sine, 


+ The upper sign always corresponds to the positive sign in (11). For flow with axial 
symmetry the x-axis is taken to be the axis of symmetry, and r = y will denote the distance 
‘rom it. Incidentally, (15) and (16) show that a Mach line along which the state of motion 
8 constant must be a straight line. 

5092.2 


P 
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this becomes 


.dr sin @si ; 
dx dp j- , ae ane along the Mach lines (r + 0), 
(17) 


where 7 denotes the distance from the axis and where j 0 for plane 


Sin pe COS 
2 r sin(@-p) 


flow, and j = 1 for three-dimensional flow with axial symmetry. 
The slope of the streamlines is 
dy/dx = tan@, (18) 
according to equation (12), and (12a) and (125) can be used as they stand: 
ridaatied along the streamlines. (19) 
wdw--(a?/p)dp = 0 J Y (20) 


VI.1. The equations (1) and (15) to (20) allow us to integrate the 
equations of motion by means of an iterative step-by- 
step process which is due to Massau (7). Suppose that 
the state of motion is known at two points, 1 and 2, 
say, which do not lie on the same Mach line or stream 
line and which are sufficiently close to one another to 
ensure that none of the variables changes greatly overa 
length of the order of their distance apart. The point of 
intersection, 3, of the Mach lines through 1 and 2 (Fig. 5) 





Fic. 5. 


may be found, to a first approximation, by replacing 
the Mach lines by their tangents at 1 and 2, the slope of which is given 
by (16). The value of 6 and p, at 3, will then be determined, by means 
of the two equations (17), from the state of motion at 1 and 2. Th 
streamline through 3 may then be approximated by its tangent, and (18) 
will permit us to determine its point of intersection, m, with the straight 
line 1-2. The state of motion at m may be found by interpolation 
between those at 1 and 2. According to (19), 8, = s,,, and p, and ag can 
be calculated from the equation of state, (1); ws will then be given by 
(20), ws by (15). 

With these first approximations for the state of motion at 3, better 
interpolated values of the variables may be used, instead of their values 
at 1, 2, and m, for a repetition of the first step, and so forth. 

A slight modification of this procedure will be necessary when the 
computation starts from a Mach line or a streamline. It must be borne 
in mind that a single Mach line (or streamline), with the state of motion 
along it, does not determine the flow. Some data must be given along 
some other line which meets it, and it is at this corner that the computa- 
tion can be started. As a typical case, let the state of motion along 4 


+ m may lie outside 1-2; the state of motion at m can then only be found if this point 
lies inside the smallest characteristic polygon corresponding to the initial data. 
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Mach line, ., of the plus-family and the shape of a streamline, S, be 
given, for example, a solid boundary, which meets the Mach line at an 
acute angle at 1, say (Fig. 6). If 2 is a point on the Mach line close to 1, 
and 3 the point of intersection of the minus-Mach line through 2 with the 
streamline, S, then the second equation (16) fixes the position of the point 





Fic. 6. 


3 and the second equation (17) determines ps, since 6, is known. 83, ps, 3, 
and w, are found as before. 

To proceed with the construction of the network note that 3-4, 5-7, 
6-8 are all ordinary lines so that the points 5, 8, 9 can be found as set 
out in the beginning of this section, whereas 6 and 10 can be found in the 
same way as 3. 

VI.2. In the special case of a perfect gas with constant specific heat, 

p CUpvexp(s/c,,), 
and we have @p/ép = a® = yp/p, ép/és = p/c,. Bernoulli’s equation, (13), 
can then be written 
H = a®/(y—1)+w?/2 = (w?/2)[14+-{2/(y—1)}sin*y] (21 
constant along a streamline, sie 
and a*p, in (17), is replaced by yp so that the density does not enter any 
more into the characteristic equations. 

In many cases of practical importance the flow starts from a state of 
uniform total energy, and H will then be constant throughout the fluid 
since it does not change even in a transition through a shock-wave. This 
is the case which was first treated by Guderley. In such a case the 
numerical process described in § VI.1 may be simplified, since 

i = a?/(y—1), T = a*/y(c,—¢,), 
da ds a 


l 2 
dp di T ds a2! — | = dQ 


p ly-L a e,—e)) 
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R. 


(for a perfect gas). This expression may be introduced into (17), and (21) 


used instead of (20), so that p and p are eliminated from the system of 


formulae used for the step-by-step process, while the new auxiliary 
variable of state, P 
Q - log, a*—8/Cy, 
is introduced into it. 
If, in addition, the flow is isentropic (and hence also irrotational, by 
virtue of (14)) the streamlines lose their characteristic properties and 
dp 


p 


wdw. 


Since H is constant everywhere, w is a function of uw only, by virtue of 


(21), and (17) can be written 


d0+-cot P wo (0-0) a, dy sin 0 nl oll (29) 
= w "7 sin(O-p) 


w 


. 
’ 


l 
where t= | cot L- <i p+kcot-(k tan p)— 3a (23) 
w 


ny 
(ls 


(cf. ref. 4), k® = (y+1)/(y—1), and a, = VH/k is the value of the speed 
of sound when w = a. The equations (16), (22), and (23) then suffice to 
determine the flow pattern: if » and @ are known at 1 and 2 (Fig. 5), the 
equations (16) determine the position of 3 as before, the equations (22) 
give f, and 65, and ps is found from (23). 

Finally, in the case of plane, irrotational, isentropic flow of a perfect 
gas (Busemann’s original case), 7 = 0 and equation (22) can be integrated 


- 6--t(u) = constant. (24) 


The velocity components are connected along the Mach lines by a relation 
which does not depend on the particular flow pattern but only on the total 
energy, H, and the ratio of the specific heats, y. The curves in the hodo- 
graph plane which correspond to the Mach lines are epicycloids (for a 
proof ef., for example, ref. 4). Moreover, since the angle between the 
velocity direction and the tangent to a hodograph curve w = w(6) is 
d = tan—(w dé/dw) 

(Fig. 7), the tangents to the hodograph curves which correspond to the 
Mach lines and the direction of the x-axis include the angle 


6+¢4 = 6+tan-( cot nv) = 6+90°+p, 
whereas the Mach lines at the corresponding point make angles 


e = OFp 
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(2!) | with the a-axis. Hence, at corresponding points, the plus-Mach line is 

mM Ol | perpendicular to the hodograph curve corresponding to the minus-Mach 
liary line, and vice versa. 

For these reasons the ‘hodograph characteristics’ have been extensively 

ised for the graphical solution of plane, irrotational, isentropic flow 

problems (Busemann’s method). It should be noted, however, that for 











» by practical purposes Temple’s alternative method (ref. 4) is more accurate 
™ 
99 
99 | 
- Fic. 7 
speed nd involves less work. It is a special case of the Massau process, the 
are ain simplification being introduced by the fact that the hodograph 
. the haracteristics do not depend on the particular flow problem. By virtue 
(92) | of (23) and (24), @ is a function of » only along any particular Mach line, 
nd the individual Mach lines of either family are distinguished by the 
ale value of the constant of integration. At any point, therefore, @ and jp, and 
voted hence the directions of the Mach lines, depend only on the ‘names’ of the 
; Mach lines which meet there, but not on the coordinates of the point, as, 
a forexample, in axially-symmetrical flow. Once the relationship between 
ation | Gandy has been tabulated, the mean slope of the Mach lines 1-3 and 2-3 
total | Fig. 5) can be calculated immediately from the state of motion at 1 and 2, 
hodo- | and the whole process then reduces to straightforward use of equation 
(fora {| (16) and the table without iteration of the steps. 
n th ; 
a VII. Shock-waves 
| Since the equations governing the change of velocity and state during 
the transition through an oblique shock-wave are well known, the method 
o the { characteristics as explained above enables us also to treat problems 


flow involving stationary shock-waves with supersonic flow down- 
stream. 

lo show howy this can be done let us consider the simple case— 

[he form in which the method of numerical integration is described in this section 


not the most profitable one for computational purposes but has been chosen rather with 


view to simplifying the understanding of the steps involved. 
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equivalent to Guderley’s original example—of a plane steady flow along 
a wall with a corner which implies a discontinuous compression (Fig. 8),+ 
For simplicity, let the gas be perfect, the wall before the corner straight, 


and the oncoming flow uniform (wy) = const., 6) = 0, H = const., s, = 0), 





At the point 1 the velocity changes abruptly (6, > 0) and hence a shock- 
wave starts from this point, provided @, is not too big. Near the corner, 
the wall and the shock-wave may be approximated by their respective 
tangents so that we have the case of an oblique shock front with known 
state of motion upstream and known velocity direction immediately 
downstream but unknown slope of the shock-wave. 

From the shock-wave equations (ref. 13) we may find the magnitude, 
w,, of the velocity immediately below the shock, the angle of inclination, 
x,, of the shock-wave, and the pressure and density just below the wave 
at the point 1. 

If the shock front can be approximated by its initial tangent up to a 
point 2, the state of motion at this point will be the same as that at 1. 
Hence, if w, > a,, there will be two Mach lines through 2, on the down- 
stream side of the shock. Taking the one which meets the wall at a point 
lower downstream we find that Massau’s process allows us to determine 
the state of motion at this point, 3. As long as the shock is extrapolated 
with constant obliquity the shock-wave equations give a constant state 
of motion just below it. If this was the case up to the point of inter- 
section, 4, between the second characteristic through 3 and the shock 
front, and if the wall was not straight from 1 to 3, so that 6, 4 @,, then 


It will be noted that according to the uniqueness theorem the flow ought to be uniform 
down to the Mach line which slopes obliquely downstream through the point 1. It must be 


»membered, however, that we were considering only flow patterns with discontinuities 
of the derivatives of the velocity components and the variables of state, in the previous 
sections. In fact, it can be deduced from the shock-wave equations that finite dis- 
continuities of the pressure travel at a higher speed, with respect to the fluid, than the 


S} yeed of sound. 
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the state of motion at 4 according to the shock-wave equations would, 
in general, not fit the characteristic equation (17) for the change in the 
state of motion from 3 to 4. 

In order to extrapolate the shock-wave correctly beyond 2 we may, for 
example, proceed as follows (ref. 7). Estimating 6, we shall get wy, a4, 94, p4 
from the shock-wave equations, and a, from (21). Equation (17) will read 


6,—0,+-0(p4—Ps3) 0 (25) 


if, for example, the Mach line 3-4 belongs to the minus-family, o denoting 
the mean value of (sin « cos 1)/(a*p) between 3 and 4. This equation yields 
, new value of 6, which may serve for an iteration of the step. In this 
way, since the oncoming flow was assumed to be uniform, the inclination 
of the shock-wave and the state of motion immediately below it, at 4, 
nay be computed. The coordinates of the point 4 are not involved in the 
formulae used and need not be determined until after the last iteration. 
In this way, we may extrapolate the shock-wave farther and farther, 
t the same time the Mach-line network downstream. 


( mstructing F 

If the oncoming flow is not uniform the state of motion at 4 depends 
on the coordinates of this point and it is necessary to resort to a more 
cumbersome trial-and-error process. We may, for example, begin by 
estimating a,. Then x, and y, can be calculated from the mean slope of 
the shock front between 2 and 4 and the slope (at 3) of the Mach line 3-4. 
The state of motion at 4, immediately ahead of the shock-wave, can now 
be found from the Mach network of the oncoming flow, and the state of 
motion just below it can be determined from the shock-wave equations. 
If the equation (25) for the Mach line 3-4 is not satisfied by the values 
thus arrived at, we shall have to try again with a new estimate for «,. If it 
is satisfied, it may still be desirable to compute x, and y,, using the mean 
slope of the Mach line 3-4, and to repeat the process with a new estimate 
for x, until all the equations are satisfied. 

For problems of axially-symmetrical flow the procedure is similar. 
Except in the neighbourhood of the axis, the radii of curvature of the 
shock-wave will be very big as against its thickness in most cases of 
practical importance, and one will be justified in computing the transition 
through the shock-wave by means of the equations for a two-dimensional 
blique shock. However, since the radius enters explicitly into equation 

7), a trial-and-error method like that described in the preceding para- 
graph has to be applied even if the oncoming flow is uniform. In the 
ase where point 1 (Fig. 8) is the tip of a body of revolution with attached 
shock-wave, Taylor and Maccoll’s solution (14) for the flow round a cone 
may be used to determine the initial slope of the shock front and the state 
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of motion at 1 and 2. If the curvature of the meridian of the body near 
the tip is not too big, one may assume that the part of the shock-wave 


which differs appreciably from the initial cone is already far enough away 
from the axis for the thickness of the shock-wave to be small against its 


radii of curvature. 


VIII. Unsteady flow 

Problems of compressible flow depending on time and one more inde- 
pendent variable can be treated in a way similar to that applied to steady 
supersonic flow in §III. There are three cases: one-dimensional unsteady 
flow, unsteady flow with spherical symmetry, and unsteady flow with 
cylindrical symmetry, the state of motion depending only on the distance. 
«, from the centre of symmetry and the axis, respectively (ref. 5). 

[f wu denotes the velocity in the direction of increasing x, and p the 
density, the equation of continuity is 


c cu Cp ae 
ea 1 Ff Ip—; (26) 
ct Ox Cx “= 

where 7 = 0 for one-dimensional flow, 7 = 1 for flow with cylindrical 


symmetry, and j = 2 for flow with spherical symmetry. If viscosity and 
heat conduction are neglected, the flow will be isentropic (except for the 
transition through a shock-wave), and by virtue of the equation of state, 
(1), the equation of momentum is 
eu cu. a* dp 2 
+ u ak F 0. (27) 
ct Or p Ox 
Introducing the new variable of state, 
is 
a 
Ww | - dp, 
/ P 
0 
the equations of motion are transformed into 
c 


+ (uta) rm ta) —ja~, (28) 
| ot ex| x 


or 


d(w+u) £(w4 u) dt +- (w+u) da 


c ox 
- {dae—(u-+a) dt} (w+-u)—ja dt, (29a) 
: “ Ox x 
d(w—u) fdx—(u- a) dt} — (w—u) ja~ dt. (29b) 
" OX x 


If w and a (and hence also w) are given as arbitrary, continuous functions 
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long some line in the (x, ¢)-plane, (29a) and (28) determine é(w+-w)/éx 
nd é(w-+-u) ct at any point of that line, unless its direction is such that 

r—(uta)dt 0. (30a) 
which case they remain undetermined and (29a) becomes 
a - au da 
l(w—+-w) lp+du ja—dt J (31a) 
{ : r u ae 


hich is a condition for the derivative along the line itself. Similarly, 
scontinuities of the normal derivatives of (w—w) in a flow in which wu 


] 


da are continuous can exist on, and only on, the lines for which 

da (u—a)dt 0: (305) 
lany given distribution of the variables, w and a, along such a line is 
mpatible with the equations of motion only if the equation 
jau dx 


lt (315) 


a, jau 
ad, . f 
p x Uu—a x 


s satisfied. 


Hence, discontinuities of the derivatives of (w-+-w) travel with the velo- 

u-+-a), discontinuities of the derivatives of (w—u) with the velocity 

a). In general, small disturbances are propagated with the velocity 
dx/dt ua, 


4 


is, with the local speed of sound with respect to the moving fluid. 
It will be noted that the Mach lines, (30a) and (305), are real curves 


rsubsonic as well as for supersonic flow. 


VIII.1. Differentiating (28) with respect to x and t, we find that 


c . . . . . 
(u-La (oo u) terms containing derivatives of lower 
ct ( i 
order, 
c oe” | bee ‘ ' : 
uta) (wtu terms containing derivatives of lower 
f 1 n| 5 
aneni ox 
order, 
oN 
d”(wtu) fda —(u--a) dt\” (w+w)- 
ca" 


that the results of the preceding paragraph hold for the derivatives of 


iy order. Hence, if the flow is uniform in one region but not in others, 
this region must be bounded by Mach lines. With much the same argu- 


ents as for steady flow, in § V, it can be shown that if w and a are given 
a segment of an ‘ordinary line’ in the (x, ¢)-plane (for example, at 
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some point, during a specified interval of time) the corresponding Mach 
quadrangle will show where, and at what time, this condition uniquely 
determines the flow. 

The integration of the characteristic equations, (30) and (31), by 
Massau’s step-by-step process is exactly analogous to that of the corre- 
sponding equations, (16) and (22), in the case of steady, irrotational, 
isentropic flow. 

VIII.2. Since the speed of sound depends on the state of the fluid only, 
the equations (31 a), (316) may be integrated, in the case of one-dimensional 
flow (j = 0), without reference to any particular flow pattern. For a 
perfect gas and isentropic flow, 

a’ = yp|/p = A*pY- 
2A 2a 
so that pry-Diz4 y a Pe u 2m constant (32) 
/ 
are the curves in the density-velocity plane which correspond to the Mach 
lines. 

The orthogonality relation between the four families of curves in the 
case of steady plane irrotational isentropic flow has its analogy for un- 
steady one-dimensional motion. If we put (ref. 5) 


-q = fu?+i = 4u?4 k- du?+ =<_2 
J p Jp 
(so that the equation of momentum, (27), reduces to 
Cq/ Cx cu) ct, 
and we may write q = e¢/ét. u = 04/02), 
the equations (31a), (316), with j 0, become 
dqt+-(u+ta)du 0. (33) 


Comparison with the equations (30a), (306) shows that, at corresponding 
points of the (x, ¢)- and (q, w)-planes, the Mach lines of one family are per- 
pendicular to the (q, w)-lines corresponding to the other family of Mach lines. 


For a perfect gas ; 
S t= a*/(y—1) q— hu’, 


and (32) ean be written 


g + wt (y—1)m( mu 0. 
| 4 } ) ) 


This equation represents two families of parabolas which correspond to 
3usemann’s epicycloids. 

This paper was written on the suggestion of Professor 8. Goldstein, and 
the author wishes to thank Professor Goldstein very sincerely for his 


criticism and encouragement. 
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APPENDIX |] 
The question suggests itself whether the systems of orthogonal coordinates 
troduced in this paper are always regular systems of coordinates, that is, whether the 
nsformation of the region of flow in the (a, y)-plane into the (a, B -plane is every- 
re a one-one transformation. Little trouble is to be expected when the orthogonal 
rdinate system is based on smooth ordinary line as one of the a-lines, for in 
t case we are only interested in the neighbourhood of this line. But it may be 


terest to find the conditions under which a regular coordinate system is obtained 


lentifying the a-lines with the whole family of the streamlines or with one of 
Mach-line families. For this purpose we have to investigate the conditions under 
the Jacobian, p 
[Po ee (34) 
o(a. y) 
fixes ind y denote partial differentiation), is bounded, and does not vanish, 


ighout the region of flow. The Jacobian is the ratio of the areas of an element 


n the plane and the corresponding element in the (7, y)-plane, and hence 

B ] hh = (35) 

a) constant is the equation of the streamlines, we can choose B in such 

way that the difference of its values on any two streamlines is equal to the mass 

w between them. Since the flow is supersonic there are no stagnation points; 

{ since the state of absolute vacuum cannot be attained in practice it follows 

t hg can neither become infinite nor vanish. We shall also exclude the possibility 

I the streamlines having corners. These are connected either with shock-waves, 

which our theory does not apply, or with Prandtl-Meyer expansions, which are 

leed singularities of the orthogonal coordinate systems introduced in this paper, 
t of a well-known kind. 

If B(x, y constant is the equation of the Mach lines of either family, hg cannot 

or be e infinite on account of two «a-lines intersecting under a finite angle 

cause the directions of the Mach lines at every point are uniquely determined 

magnitude and direction of the velocity; hg can vanish at an envelope of the 

M lines, but this appears again to be connected with the presence of a shock- 
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The only condition to which a(x, y) is subject is the orthogonality relation 
x,,B x, P 0. (2 


We are free to impose an additional condition in order to determine a. Let us choos: 


the condition that 


h, l (36 
on one of the a-lines. 

It ensures the regularity of the coordinate system on the a-line where this 
condition was stipulated, but it does not immediately exclude the possibility of ar 
envelope or a ‘thinning-out’ of the normals in some other part of the field of flow, 
However, it was shown in $ IT (3), that 

oh ,/ep h he k,,> 


and therefore h, exp | xk, hgdB, 


the integral being evaluated along a normal (a constant). Hence, if J and its 
reciprocal are bounded at one point of a normal they are bounded all along it unless 
kK, is singular at some point of the normal of an order equal to, or higher than, th 
reciprocal of the distance from that point. 

Still, there need not exist a Mach line (or streamline) which meets more than a 
part of the f-lines in the field of flow (in axially-symmetrical flow, for example, all 
the Mach lines end on the axis). In order to extend our argument to the other B-lines 
we may proceed as follows. Suppose that a m on the a-line where (36) was stipu- 
lated, so that this a-line does not meet the B-lines on which « m. Now consider 
another a-line on which a n, say, with n m, and which meets the §-line on which 
, m. The condition (36) for the original a-line does not determine a on that 
part of the new a-line where m < « n, and we are free to postulate the same 
condition for this part of ‘the new a-line. We thereby ensure the regularity of th 
coordinate system in the part of the field of flow where m ’ n, provided x 
satisfies the condition mentioned above in this region. Repeating this process we 
can extend our argument to any part of the field of flow if we assume that the 
region of flow is singly connected—otherwise we should hardly expect a shock-free 
flow pattern to be realized in practice. 

The condition for the curvature of the «a-lines is not altogether without interest, 
in spite of our assumption that the a-lines have no corners. In fact, it will be shown 
in Part II that singularities of the curvature of the a-lines ean occur on the axis, 
in axially-symmetrical flow, which do not imply the formation of a shock-wave no! 


destroy the regularity of the systems of orthogonal coordinates. 


APPENDIX II 
NOTE ON THE CALCULATION OF AXTALLY-SYMMETRICAL 
SUPERSONIC FLOWS 
By 8S. GOLDSTEIN 
(Department of Mathematics, The University, Manchester) 

It appears from Dr. Meyer’s account that, as a result of recent German work, 
together with the improvements in presentation which he has himself introduced, 
the theory of the method of characteristics for axially-symmetrical supersonic flows 
is now in a form which is both useful and easily comprehensible. One point only 
still seems to need clarification, namely, conditions at and near the axis. To start 
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ith we may put a definite question—how shall we calculate a supersonic, axially- 
nmetrical, irrotational, isentropic flow when the velocity is given along the axis ? 
th the velocity given along the axis the local speed of sound is known from 

Bernoulli’s equation, and the pressure and density from the adiabatic condition. 


We denote distance along the axis by x and perpendicular to it by r, and the 


locity components in the directions of x and r increasing by wu and v. 


The equation of continuity Is 
c Iv 
(pu) (pv) - f QO: (1) 
‘ cr r 
0 when 0, so from (1) 
v » l/cu ucp 
lim et Sar, a r} (2) 
07 CY] =o 2\OX P CX p= 
reover, the axis is a streamline, and therefore 


wdu (1 p)dp udu- (a?/p)dp 0 


the axis. since u w and ds 0. It follows that 
° 
u cp u~ cu 
(— —f) (——) , 
p Ca 0 aQa* CL r=0 
, v 1 cu . » cu P 
nd henes lim (M2—1) $ cot®?u—. (3) 
? Z ON = Cx 
Hence, at a small distance dr from the axis, 
, eu 
sor cot™ - (4) 
Pe Cx 
ypproximately 
The equati n for the absence of vorticity is 
ct cu a 
0, (5) 
CA cl 
1 since v 0 on the axis, Cu/e7 0 on the axis, and at a small distance or from 


axis u is unchanged to the first order. 
We have now to check that these results agree with the results obtained from the 
thod of characteristics. Dr. Meyer’s characteristic equation (17) is not, as it 
stands, valid at a point on the axis, but we may find its limiting form. 
Since the flow is irrotational and isentropic, dp/p wdw. Also 
a w sin p, (6) 


d the symbol R introduced by Dr. Meyer is equal to v/r. Hence the equation 
wv , - 
dw +w tan pd@—- tan psin ph, da 0. (7) 

= 


upper or lower sign is to be taken according as the velocity component, vg,, 
rmal to the Mach line is positive or negative; note also that h,da must be given 
ts correct sign in (7). 
On the axis the characteristic equation becomes 


Cu 
dw -+-w tan w.dé—}cosp h, da 0. (8) 
c 


a 
Let A and B be two points on the axis at a distance 6a apart, and take the Mach 
nes AC, CB as shown in Fig. 9. If 8 is sufficiently small, the variation in pw over 
the distance 62 may be neglected for a first approximation. Positive directions are 
licated by arrows. Denote values at A, B, C by the suffixes 1, 2, 3. 
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The distance, 6r, of C from the axis is da tan p. Hence by our previous result Dr. } 

Cu. entro 

Vs H cot pu De Ox 1 (Ws, — Ww, )cot pe (9) | Sin 

(since u won the axis). do th 

On AC, vz is negative, and h,da is positive, and equal to 362 sec, from A to C, tion { 
Hence W,—W, — Ww, tan up— }(w,—w,) = 0. (10) 


(We have again replaced (éu/éx)dx by w.—.) 








and « 
' 
Hence 
and 
M so th 
If 
Fic. 9. ; 
Pip 
On BC, vg is positive, and h, da is negative, and equal to — }éxsecp, from B to 
C. Hence 
W,—W,+w, tant }(w,—w,) = 0. (11) 
= : . 5 ; . s,m 
The solution of these equations for ws and 6, is — 
entre 
, , , . ; Lies . 9) 
W, = 3(w,+w,), w, tan, = }(w y—w)). (12 total 
Since v, = w63, and to the first order wy us, these results are consistent with | Sir 
our previous results. If we take sufficiently small steps along the axis, we may 6, ar 
therefore start the calculation from the axis in this way. For w and yw we should | only 
take the values at the mid-point of AB. f the 
[t is clear that for the initial step the method of characteristics has no advantages | equa 
over more elementary methods. If the velocity is given analytically along the axis, entr 
so that its derivatives are easily and accurately found, then away from the neigh- a sa! 
bourhood of any singularity in the given distribution along the axis, the velocity fully 
at a small distance from the axis is probably best found by assuming expansions 
of the form . ‘ 
u a,+b,r?+¢,r 
p a,+b,r?+cyrt+... : (13) 
” Ay r- bs r j Cs r? 1 
where the a, b, c are functions of x, substituting in the equation for the absence | 


of vorticity, the equation of continuity, and Bernoulli’s equation, and equating | 
coefficients of powers of 7, The functions a,, a, are assumed to be known functions 
satisfying Bernoulli’s equation, and then az, 6,, b,, ete., may be found in turn. 

If we are calculating the flow through a nozzle with a given boundary, we arriv' 
at a step where the conditions at A and C are known. The equations (12) then 
provide a check on the accuracy of the values at A and C (since we must have 
3— WU, = 2w, 6, tan (14 
nearly), and enable the value, w,, of w at B to be calculated. The equation (14) has 
been compared to the results of a numerical example where the variation of thi 
velocity on the axis was particularly strong, and it has been found well confirmed. 

In the more general case of a flow which is neither irrotational nor isentropic th 
motion is not uniquely determined by giving only the distribution of the velocity 
on the axis. We need to know such additional data as will, in conjunction with 
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7 | D: Meyer’s equation (14), determine the distribution of the total energy, the 
tropy, ana t! vorticity. 
i) Since the axis is a streamline the equations (1) to (4) hold also in this case, and so 
the equations (6) and (9). Equation (5) has to be replaced by Dr. Meyer’s equa- 
1 1 I ‘ A 1 
‘ n (14), or an equivalent equation. Instead of equation (7) we have 
1 to ( | { 1 
} > ) v ° - 
10 dp pw*tan pad pw sin p tan ph, do 0, (7a) 
lon the axis this becomes 
cu 
dp pu “tal dé Spw cos — h, da 0, (Sa) 
Cr 
nee P3—P,1+p *tan 6, [pw(ws, w,) 0 (10a) 
Ps Po pw*tan 6, [pw(w, w,) 0, (lla) 
Pp Pe 
at 4(p1+ Po), w tan pO, > —1(w,—w)). (12a) 
2pw 
f we take for w, p, x the values at the mid-point of AB and if we remember that 
vdw 0 on the axis, which is a streamline, we come back to 
nB 1 U3 1 > 
LO tan pL t\We2 Y1)> Ws w o(W, T U's). (12) 


6, 
ust be determined from the additional information about the distribution of the 


tropy, and p, and a, can then be found from the equation of state, provided the 


tal energy at C is known. 

t wit Since 6w/ér, Op/ér, 6u/er do not vanish on the axis, the approximation (12) for 
ve mi §, and ws, is not quite as good as in the case of irrotational, isentropic flow where 
shot nly terms of the order of 62? are neglected in these equations. In the general case 
neglected terms are of the order of (€u/ér)da. This applies also to the ‘check’ 

untages juation (14). However, in most practical cases it is found that the variation of the 
he axis ntropy is very small and then we shall expect the equations (12) and (14) to provide 
neigh | satisfactory approximation also in the more general case where the flow is not 


elocit fully irrotational and isentropi 
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By E. T. GOODWIN and J. STATON 
(The National Physical Laboratory, Teddington, Middlesex) 
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SUMMARY 
This table gives five decimal values of J(jo,,7) for nm = 1(1)10 and the rang 


r 0-00(0-01)1-00, 75, being the nth zero of the Bessel function J). It is useful 


in the solution of some electrostatic problems with axial symmetry. 


A SINGLE-VALUED solution of Laplace’s equation in. cylindrical polar 
coordinates (r,¢,z) is formed by linear combination of terms of the form 
J (krjexp(+kze+ind) and Y, (krjexp(+kz+tind). If the solution is to be 
symmetrical about the axis r = 0, then only zero-order Bessel functions 
will appear, and if it is also required to be finite on the axis, then the 
linear combination will consist of terms of the form Jj(kr)exp(+kz) only. 
Suppose now that the solution V is to take a constant value V, over the 
surface of the cylinder r = 1, then at points inside the cylinder 
=} » 

V = H+ 214 n EXP( Jo,» z)+ B, exp( —Jen 2) Wo(jon r), 
where jp, is the mth zero of the Bessel function J. 

In electrostatic and other problems with axial symmetry the solution 
may frequently be obtained by fitting expansions of the type given above 
at the appropriate boundaries. Professor D. R. Hartree, who was at the 
time interested in just such a problem, suggested to us that the preparation 
of a table of Jo(jo,,7), for values of r in the range 0 to 1 and the first few 
values of n, would greatly facilitate such work. 

In the accompanying table five decimal values of Jo(jo,,, 7) are given for 
n 1(1)10 and the range r = 0-00(0-01)1-00. These values were calcu- 


] 


lated to seven decimal places by interpolation from the B.A. Mathematica 
Tables, vol. vi, of Bessel functions; they were then differenced in the 
r direction on the National Accounting Machine as a check and all 
doubtful roundings off were examined; all five decimals printed should, 
accordingly, be accurate. 

The work described above has been carried out as part of the research 
programme of the National Physical Laboratory, and this paper is 
published by permission of the Director of the Laboratory. 
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THE APPROXIMATE EXTERNAL AND INTERNAL 
FLOW PAST A QUASI-CYLINDRICAL TUBE 
MOVING AT SUPERSONIC SPEEDS 
By G. N. WARD 
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SUMMARY 

By using the Heaviside operational method applied to the linear perturbation 
eory for supersonic flow past a cylindrical tube of nearly constant radius at small 

neidence, the aerodynamic force and moment coefficients resulting from the fluid 
pressure have been calculated for the outside of the tube, and inside the tube for 
a length (1/2—1)* times the diameter of the tube, where M, is the Mach number of 
the undisturbed flow. The results are obtained in a form suitable for computation. 
The expressions for the drag, lift, and moment coefficients involve only the slope 
f the meridian section at army point, the angle of incidence, and certain auxiliary 
functions which are the same for all tube shapes. 

The distribution of the velocity components inside the tube has been considered 
nlinear theory ; singularities and discontinuities occur in these velocity components, 
r in their derivatives, corresponding to a discontinuous change in the slope or 
urvature respectively of the meridian section of the tube, or, at entry, of the 
streamlines bounding the fluid entering the tube. The nature of these singularities 
is studied. 


| 1. Introduction and mathematical formulation 

Tue Heaviside operational method (1) will be used to determine the 
approximate external and internal supersonic flow past a tube of nearly 
constant radius at small incidence. Lighthill (2) has obtained the formula 
for the external drag, and the present author (3) the formula for the 
external lift by using more laborious methods. Certain auxiliary functions, 
which occur in the expressions for the pressure, also occur in the theory 
of bodies of revolution whose meridian sections have discontinuities in 
slope (4). 

Consider supersonic flow past a tube of nearly constant radius R whose 
section has small slope and whose axis is at a small angle « to the stream. 
Take a system of cylindrical polar coordinates, r,, 0, x such that the x-axis 
is along the tube axis and x = 0 at the mouth of the tube, and let &(2) be 
the slope of the tube section at a distance x from the mouth. Let U be the 
velocity of the undisturbed stream, which will then have components 
Ucosa parallel to the tube axis and U sina normal to the tube axis; we 
choose @ = 0 to be in the direction of the component U sina. If v,, vg, v, 


are the components of fluid velocity at any point in the directions of r,, 0, 
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and 2 increasing respectively, then the boundary condition of zero normal 
velocity relative to the tube is satisfied if 
v, = §(x)v, (1) 
at the surface of the tube. Since the flow is assumed to be supersonic, 
disturbances cannot be propagated upstream, and we must have 
zero disturbance velocity (ie. v, = Usinacos@, vg = —Usinasin§, 
v, = Ucosa) when x < 0. 
We assume that the fluid is a perfect gas and that the flow is isentropic 
and irrotational; we shall neglect viscosity and heat conduction, so that 
the disturbance velocities can be derived from a scalar potential U®, 


and we may write pace ea 
; a. , = U(sinacos 0+ 0e/ér,), 


vg = U(—sinasin 6+ 6®/r, 26), (2 


v, = U(cosa+éO/éz). 
Since the slope of the meridian section of the tube and the incidence are 
small, the disturbance velocities will be small compared with U and we 
may neglect the squares and higher powers of the disturbance velocity 
components in the equation of motion. In this way we obtain the well- 
known Prandtl—Glauert equation for ® 
emp led 1eo oD , . . 
ot, ot oe age = MOH) " 
1 41% 1%} 
where M, is the Mach number for the undisturbed stream, and the boundary 
condition (1) becomes 
U(acos 0+ 6®/dr,) = Us(x). (4) 
Equation (4) may be split into two equations, 
o®, /or, = 3(x), o®,/er, = —acos8, (5) 
where ® = ®,+9,; ®, is the potential for zero incidence and Q, is the 
extra potential due to incidence. Since the radius of the tube is nearly 
constant, the boundary condition (4) or (5) may be applied on the cylinder 


r, = R, x >0 without causing loss of accuracy. If AP is the extra 

pressure due to the disturbance, then Bernoulli’s equation reduces to 
AP oD 
— = ——, (6) 
pU? Ox 


where p is the density of the undisturbed stream. 

In order to simplify the working slightly, we make the transformations 
«= BRz, r, = Rr, 0= 8 and write O = Rd, 3(x) = 3(BRz) = (2), 80 
that equation (3) becomes 

Ab Led 1d _ a Q) 


or2 " rér ° 72062 = az? 








and 
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sin 
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col 


si 
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and the boundary conditions (5) become 








ormal 
66/ér = s(z2), 0¢/er = —acos8, (8) 
(] when r= 1,2 >0. 
sonic. The disturbance pressure is now given by 
have AP 1 o¢ (9) 
sin 8 pU2— =Béz- 
If p is the Heaviside operator for z, then the operational form of (7) is 
tropi 
*h 1 l & n . 
» that a 2 Oe < we) $ = p*d (10) 
| Uo. ore ror £0 
since Ob /ez 0 for z < 0. 
The solution of (10) which is independent of @ is 
i ¢ = Al,(pr)+CK,(pr), (11) 
| where J, and K, are modified Bessel functions. 
2 are _a ° . ° ¥ 
If S(p) is the operational representation of s(z), then the boundary 
1d we a ° ° . 
oc conditions (8) become, in operational form, 
oclty 
edi. Od /ér S(p), 0d /éer —acos@ when r= 1. (12) 
From the second of (12) we see that the solution required for the extra 
potential due to the incidence will contain a factor cos@. The solution 
" of (10) satisfying this requirement is 
dary d = |A, ](pr)+C, K,(pr)|cos 6. (13) 
The equation (9) for the pressure becomes 
(4 a 
AP/pU? = —pd/B, (14) 
, since d = 0 for z = 0. 
vo 
sthe | 2, Solution for the external flow at zero incidence 
early As may be seen from the asymptotic expansions for the Bessel functions, 
inde1 the term J, in (11) represents waves travelling outwards upstream, so that, 
xtra since ¢ and its derivatives must vanish for z < 0, we must have A = 0. 
° K, represents waves travelling outwards downstream and so is appropriate 
P for this case. 
From (11) and (12) we have 
pCK9(p) pCK,(p) = S(p), (15) 
b10ns 
K,( pr 
), 80 and so o= — sol Pr) S(p); (16) 
pk,(p) 
( f K r 7 
- also r = — ol Pr) S(p). (17) 


; Oz K,(p) 
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On the tube, this gives 





AP 1 K,(p) g 
= — S(p). 
U2 B K,(p) (P) (18) 


p 
The pressure for small z (> 0) is 








AP(0) . 1 KA) ¢ s(0) s(0) 

— = lm = S— SQ) = —- = 

pU® ~~ BE) ~ B ~ (el) (19) 
in agreement with the two-dimensional result due to Ackeret. 

We shall now use the product theorem: if 
filz) = Kp), f2(z) = Ep), 
then alee) ues | AMAe—0 dt. (20) 
P 
0 


This must be used to interpret (18), since if F(A) is put equal to AK, (A)/K,(A) 
in the integral one 
aif 


fz. =— [ Fa) d 1K, (21) 


2710 





J 


L 
where ¥ is the contour c—ioo to c+-ico passing on the right of all the poles 
of the integrand, this integral does not converge. 

Consideration of (19) leads us to write (18) in the form 


AP Pe | K,(p)—Ky(p)) 
- — _¥§ i, ape oR, 22) 
pU? B (P)) K,(p) ,’ \ 
and so we have 
AP ] 4 y »)9 
pU? = 7" | WW (z—t)s(E) a}, (20) 
0 
where W (z) = PK, (p)—Ko(p)] (24) 
K,(p) 


The result (23), with W(z) defined by an integral equation, was given by 
Lighthill (2). The function W(z) was tabulated to five decimal places for 
z = 0(0-2)10 by the Admiralty Computing Service (6); the values for 
1-8(0°2)0 were computed by the present author and Miss Roussak at 
Manchester. W(z) is given in Table 1 on p. 231. The properties of W(z) 
will be considered in the next section. 
The external drag coefficient for a tube of length 1, = | BR is (cf. ref. 2) 
1 11 


ext. drag — Jf ? 
Cpovext) = jpU = R? = "7 (2) dz—2 | [n (|z—t|)s(z)s(t) dzdt. 


0 00 : 
(25) 


The velocity components away from the tube may be found from (16). 
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In particular the shape of the streamlines is given from 


op _ K,(pr) , 
= S 26 
or Kip) (p). (26) 
éd/ér gives the (approximate) slope of the streamlines at any point. 
By using the asymptotic expansion for K,, we find that @¢/ér = 0 for 
-r+] 0, and for z—r+1>0 





os 
oP we 14 O2—r + 1)] 


cor Vr 


= 0), /(F\p+oe—r+ 0} 
ry 


Thus the deflexion through the initial shock-wave decreases with increas- 


(27) 


ing r or 7,, and hence we may expect that the shock-wave will be curved 


in a more exact representation of the flow. 


3. The function W(z) 
From (24) on interpretation by (21) we have 
L [ K,(A)— Kg(A) 
) 


edz dH. 28 
Kin e” ¢ (28) 


Now the asymptotic expansions of Ky, K, are 

- (0, m) 

+ 2, (2A) m1? 

‘ 7 \4 — (1, m 

K,(A) ~ | -eAl1 4 >! 
m=1 


(n+m—4)! 


n! (n—m—}4)!" 


K,(A) ~ (=) ‘ 





(29) 


where (7, Mm) = 


Also the K’s have branch points at A = 0, and_K, has no zerosin |argA| < z. 
Hence we may transform ¥ into the contour .@ (—o, 0+). 

Let the contour 4 become the upper and lower sides of the cut in the 
\-plane along the negative real axis plus a small circle enclosing A = 0; 
then, integrating along the three parts of the contour and using the results 


K,(Ae**”) K,(A) ial, (A), 
K, (Ae*') kK, (A) ial, (A), (30) 
Ty(A).K, (A) +4 (A) Ko (A) 1/A, 

e e-uz du 


we have W(z) — eee 3] 
K?(u)+relu) u aa 
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The function W(z) is defined by (28) only for real z > 0, but we may now 
take z to be a complex variable, in which case (31) provides the continua- 
tion into the strip —2 < Rl(z) < 0 since the integral converges in this 
region, as may be seen by using the asymptotic expansions of J, and K,. 

By substituting the asymptotic expansions (29) in the integral (28) 
expanding the integrand in powers of 1/A, and integrating, we have the 
Taylor expansion wi 


2) = 4$—#2+- #2"... (32) 

which will converge for |z| < 2. 

From the integral (3 = we can determine the nature of the singularity 
+6, 


at z= —2. Putz = 














then W(—2+¢) = [ g(uje-“S du, (33) 
0 
ebu 
7} ITA ‘ 
where g(u) = Wf KX (u) tea)" (34) 
When w is small, g(u) = u+O(uPlog u). (35) 
, a. 1 a 
When ~ is large, g(u) ~ —{14+—+0|[—}}. (36) 
F | 4u u* 
Therefore 
W(z) = wu) du+ {| 1+. 5 +0(7 ‘jes du 
Z 4u u? 
0 
2 e-$ r 
= O(1)+- —+— Ei(¢)+ o( | (37) 
al Qs 
1 
9 ‘ 
= oF SF tee t4001) 
aC 27 
that is, as z > —2 
: 2 | 3 . 
W(z) = = 5 +5 log 5+ OU ). (38) 


By substituting the expansions of the Bessel functions in ascending powers 
of A into the contour integral (28) taken for the contour -%, or into the 
infinite integral (31), in which the most important contribution is from 
small u for large z, we may obtain the asymptotic expansion 
W(z) ~ ed ; (log 22—§)+ (39) 
z 4)+... 


es ' 


we w4 


for large z. 
When z = 10, this expression gives W(10) = 0-01099. 
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y now TABLE 1 
tinua- —- - — 
n this a , we wee eee. ee = 
d K. 20 x 20 O14542 | 60 0°02936 
a 18 4°896 22 “13131 | 62 | 02762 
(28 16 2825 2°4 ‘11889 = || 6-4 02602 
re the 1-4 i 2°6 ‘10798 || «6-6 "02455 
1-2 1°52291 28 09833 | 68 "0232 
(32 | 
1-0 I 507 30 0°05976 70 0°02194 
0:8 98461 372 “08216 72 "02078 
larity 06 ‘81041 3°4 "07536 74 ‘01970 
: 0-4 68644 36 “06929 | 7:6 “01869 
-O°2 *58324 38 -06384 | 78 | °01775 
aaa re) 0*500 40 0°05896 8-0 0°01688 
hata +0°2 4319¢ 4:2 05455 «|| 8-2 01608 
0-4 37552 44 705058 = || «84 01534 
0-6 *32839 46 04697 || 86 | 01463 
(34) 08 -28867 48 *04370 88 | 01397 
fe 1-0 0°25497 5:0 004073 9:0 0°01335 
vv 12 22617 5:2 03803 92 01277 
I*4 "20143 5°4 "03556 9°4 “01222 
(36) 16 *18009 56 03331 9°6 ‘O1I70 
18 16155 58 "03125 9°8 ‘OII22 
2°0 14542 6-0 0°02936 || 10-0 0°01077 
4. The result for the external flow at incidence 
The extra potential at incidence has the form given by (13). Since the 
disturbance vanishes for z < 0, we must have A, = 0. The second con- 
ate dition in (8) then gives 
pC, Ky(p) —«a. (40) 
K,( pr 
Hence dh x COS g Dal (41) 
pky(p) 
f K , 
CO ” 
and = —acos6- iPr) ‘ (42) 
(38 Oz Ky(p) 
— The extra pressure on the tube due to incidence is 
to the 


AP x.cos@ K,(p) 
TO } U = = W'7..\° 43) 
trom pU? B Kip) 


For small values of z ( 0), the pressure is 


(39 
Af 9) = lim «cos g Ail) ae a. ? (44) 
ail ee Ki) (M3—1) 


again in agreement with two-dimensional results. 
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The equation (43) may be written in the form 








AP __ «acosé V(2) , 
U2 iene Zz), (45) 
= K,(p) 
There J (z a se. 
where ) K'(p) (46) 


The properties of V(z) will be investigated in the next section. 
The tube being of length 1, = /BR, the normal force coefficient is 


l 
, normal force, (7... _— = 
Cyiext) —_ = U27R2 = 4 J (z) dz = 2aL(L), say, (47) 

0 
the normal force being in the direction of the component U sin a of the 
undisturbed stream, and the moment coefficient about the origin (that is, 
the centre of the mouth of the tube) tending to decrease the incidence is 


1 


= 2a0B | zV(z) dz = 2aBM (I), say. (48) 


« 


0 


a moment 
M(ext) = bp U 27 R3 
There is no extra drag measured in the axial direction. Further con- 
sideration of the forces will be deferred until we have investigated the 
properties of V(z). 


5. The functions V(z), L(z), and M(z) 
Interpreting the operational form of V(z) given by (46) we have 


i K,(A) ,y.dA 
277 . K4(A) yy 


£ 


K;, has an asymptotic expansion 


a a \i 1 ~— (0, m)+(2, m) 7 
BAe — 2e-AlL 4 on a 50 
— (5) | "2 pa (2a) it 


(that for K, having been given in eqn. (29)). 

By a similar process to that used for W(z), the contour Y may be trans- 
formed into the contour -%, except that in this case K}(A) has two zeros 
in a7 < |argA| <7. These are at A —0-6435-+0-50127 (6). Thus we 
obtain the result 


V(z) = e-964352{ ] -2120 cos(0-501 22) +-0-1898 sin(0-5012z)|— 


Jf iA) n.d 
2mi J K4(A) A’ 


AM 
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| and using formulae analogous to (31), we obtain finally 


V(z) = e-6435-[ 1-2120 cos(0-5012z)-+-0-1898 sin(0-5012z)]— 


(45 
[ e-u2 du (52) 
aT eee TE he v2 
(46 J KYP(u)+77142(u) u? 
=? 0 
This result is due to the Admiralty Computing Service (6). 
Again by considering z to be a complex variable, it will be seen that we 
have continued V(z) into the range —2 < RI(z) < 0, since the integral in 
52) converges for R1(z) 9. 
(47) | Expressions for V(z) analogous to those obtained for W(z) may be 
| obtained by the same methods. 
f the As z > —2, we have, 
at is 2 l 
biste V(z) = ——log 53 FOU). (53) 
| The series expansion in powers of z is 
18 
V(z) = 1—4e—}2? + ..., (54) 
which converges for |z| < 2, and the asymptotic expansion for large z is 
con . 
1 the m1. 2 , 24 a 3 ~~ 
J (Zz) a om 33 35 log ae +.... (55) 
The exponential terms in (51) or (52) are negligible theoretically, but must 
in fact be taken into account when V(z) is computed for z ~ 10, say. 
| When z = 10, this expression gives —0-00191; by including the expo- 
49 nential terms we have V(10) ~0:00163. 
(4! , ‘ ; 
For the function L(z) we have 
} T K,(p _ 
L(z) = V(z) dz = —— (P) (56) 
/ pKy(p) 
(m0 0 
[he series expansion in powers of z, most easily obtained by integrating 
54), is 
y »—_ Jo2__ 151 " 
L(z) = z—}22—H2 +... (57) 
rans- \ The asymptotic expansion for large z is 
ZeTOS 
we l Ko 
S L(z) ~ 14+—+.... (58) 


In order to obtain the operational representation _f M(z) we require the 
known result that if f(z) = F(p), then 


(51 ' d {[F(p) 
fe) = —P(} 
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Thus for M(z) we have 


= oe F a ae d K,(p) . 
M(z) = [=z dz = (ee (59) 
0 


Either from (59), or more simply from (54), we can obtain the series 


expansion , 7 
. M(z) = 42?—328— fet ..., (60) 


6 
and from (59) we get the asymptotic value for large z 
2 
M(z) ~ —-=. (61) 
A table of V(z) was computed by the Admiralty Computing Service to 
five decimal places for z = 0(0-2)10 and is given in Table 2. It will be seen 
that V(z) has a zero near z = 3: it is believed, but has not been proved, 
that this zero is at z = 7. 


TABLE 2 ‘ 








2 V (z) V(z) 2 | V(z) 
oO I*00000 40 | 04170 || 8-0 |  —o-00861 
o-2 0°89827 ] 42 | "04490 82 | “00731 
0°4 79576 || 44 | 04650 || 84 “00618 
06 | ‘69560 || 46 | 04680 } 8-6 "00521 
os | “60006 | 48 | 04607 || 88 00438 
| ‘ i] 
1] | 1} 
0 | +0°51064 | 5:0 —0'04456 || 9:0 | —0:00366 
I2 "42029 «|| ~52 *O4247 \ 9:2 "00304 
4 "35356 | a4 | 03995 || 94 | 00253 
| 4 | | | 
16 | -28661 || 5°6 03721 || 96 | ‘00211 
18 *22738 58 03431 «|| 98 | 00176 
| | 
20 | +0°17559 || 60 —0°03132 | 10°0 | —0°001 47 
} | 
22 | "13089 || 162 | 2837 | 
24 «| 709275 || 64 | 02551 ! | 
2°6 06067 6°6 *02277 | 
28 | "03405 68 | ‘o2018 || | 
| | | 
| \ 
30 + 0°01 232 7O | —O°OI775 | | 
32 «| 000507 || 7:2 | “01556 | 
34 | 01867 74 01353 || | 
36 | ‘02899 || 76 ‘o1169_— || 
38 | 03653 «|| +78 | “01006 
} | 
40 | —oo4170 || 80 | —o-0086r_ || 


/ 1] 





6. The external force coefficients 

For convenience we summarize the results obtained in the foregoing 
sections. The force coefficients usually required are the drag coefficient Cp, 
obtained from the drag force measured in the direction of the main stream, 
and the lift coefficient C,, obtained from the lift force measured in a 
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direction normal to that of the main stream. The lift and drag forces 
may be obtained from the normal force and the drag force at zero incidence 
byresolution. In this way we obtain for a tube of nearly constant radius R 


and length h D(ext f ‘DO(ext) + 207 L (1, BR), (62) 


Ciext 2aL(l,/BR), (63) 
where Cpovext) 18 given from (25). 
The moment coefficient is 
Cwrext 20 BM (l,/BR). (64) 
In many applications / = /,/BR is quite large (say, about 10 or greater), 
so that we may use the asymptotic expressions for L(z) and M(z). We 


have then 


Cnext) ™ Cporexty + 20°, (65) 
( Liext) ™ 2%; (66) 
and for the centre of pressure 
2 B? R? - 
top ~~ —-—— (67) 
L, 
for large 1,. 
7. The results for the internal flow at zero incidence 
The appropriate solution of (10) in this case is 
d Al,( pr) (68) 
ind the boundary condition (12) gives 
pAly(p) = pAl(p) = S(p). (69) 
I pr) 
Hence d = ao S(p) (70) 
pl,(p) 
Od 2 
and a LAP") gp), (71) 
Oz Lp) 
The pressure on the inside of the tube is then given by 
> 
ae ns ~~ in. (72) 
pl - B I(p) 
The pressure for small z (> 0) is 
(0 s(0 
ae °) — — lim 4 LA) sy) = — m2 = (73) 
pl « loo B D(A) /(Mj—1) 


If we proceed as before and apply the product theorem to (72), it is 
‘ound that the formula (21) gives a divergent integral for the interpreta- 
(Pp) L(P) gi 


i\P)) 


tion of p). This is due to the more complicated behaviour 
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, 
of the I’s at infinity ; for example, the asymptotic expansion for J,,(A) for | = Py 














—thn <argA < gis rigo! 
>A se . »—AH(n+4)7i = (n, m) 
ra Sig ah. — m (n ™m) alt 1 . 
nA) (27d)! | 2, ane |+ aan [+ > (2ay™ 
(74) 
We can overcome the difficulty by defining a new function s,(z) such that 
in t 
s(z) = s(0)+8,(z), 8,(z) = S,(p). (75) asvi 
Equation (72) may then be written F 
AP ie L(p)—L(p) L(p)—L(p) _¢ oe 
— == —— | S(p)— a || ie i ry}. 6 
pU? 5| f I(p) pl,(p) Poirt (79) 
Define 7'(z) by 4(p)—4olP) = T(z). (77) | and 
I(p) | inci 
All the operational representations in (76) have a meaning, and we find 
R ; 
AP F , 
=F | a(e)—s(0)7(@) — | sre» at). (78) Psat 
pvt B| 
0 | 
From (77) we have 7'(0) = 0, and hence, integrating (78) by parts, the 
| of! 
Af ‘J . F ; sing 
—- = — —{ &X(z)—8(0 z)— 8 2— 
a | s2)—a(0)T(2)— | s(0T"e—1 a| 
0 
, | pra 
a 1 ve | ee ot 79 for 
- tore | T’ (z—t)s(t) a, i) pe 
, ger 
By interpreting (77) we have SD oe 
} 
T(z) = as [ =e. (80) 8. 
Qni) LAA) r 


. . . ° e an 
If we now substitute the asymptotic expansions for J, and J,, we may : 


expand the integrand formally in powers of e~ and interpret term by term. 
The convergence or otherwise of the resulting series does not concern us; 
each term in the integrand has a factor e*-2”), say, and the corresponding wh 
term in the integral vanishes for z < 2n. Thus there are always a finite 
number of terms and the question of convergence does not arise. In 
particular, for z < 2, we may use the first series only in (74). By com- 
paring this with (29) we see that the series for 7"(z) near the origin is the of 
same as that for —W(—z). Hence we have 


T'(z) = —W(—z) (0 <z <2). (81) 
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Professor Goldstein has shown that this last result may be proved 
rigorously by substituting 


Hwy (ze!7*) | 


oz) = }H(ze-¥7*) 11 — ) | 39 

ol ) D) ( ) H®) (ze bart) ( ) 
’ (1)(sph7i 

I,(z) = =H (ze ini) fy — Meat A (83) 


2 | H(z bt)| 
in the integral (80), expanding 1/[,(z), and considering carefully the 
asymptotic behaviour of the Hankel functions. 

From (79) we have 
AP 1 | ' ) 


: s(z) 4 W (t—z)s(t) dt (oO <2 <2), (84) 
U2 BI a 
0 
nd there results an expression for the internal drag coefficient at zero 
neidence similar to that obtained for the outside, 


1 . 3 
Canines 4 | s*(z)dz+ 2] | W(—|\t—z|)s(z)s(t) dtdz, (85) 
0 00 
valid for 0 l Ze 
If T(z) is required for z 2, then (80) must be expanded by the residue 


theorem and the resulting trigonometric series computed. The behaviour 

{ T(z) is rather complicated since 7'(z) has certain periodically distributed 
sincularities and discontinuities whose forms are discussed later. 

That the result (85) is only valid for 0 </ < 2 is not a very serious 
practical disadvantage since the flow does not usually remain supersonic 
or very far down the inside of the tube, but breaks down to subsonic 
through a system of shock-waves. Moreover, it will be found that, in 
general, the pressure on the tube becomes logarithmically infinite at z = 2 


ind the approximations of the linear theory are no longer valid. 


8. The results for the internal flow at incidence 
The work is very similar to that for the external flow and we can obtain 
analogous results. 
The appropriate extra potential is 
bd = A, 1(pr)cos 8, (86) 
vhich gives, on the inside of the tube, 
AP vcos 6 L,(p) 
pU? B I'(p) 
As before, by expanding the asymptotic expansion of J,/J, in powers 
ole~’, and comparing the first term with (43), we find 
AP xr cos 6 ,, 


= V(—z <= 2 < Zp. 88 
> a ) (0<2 <2) (88) 
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and hence the normal force coefficient is 
I 
Cyant) = 2a | V(—z)dz (0 <I <2). (89) 
0 
As for the external flow, the extra potential due to incidence gives no 
axial force component. 
We thus obtain for the force coefficients parallel and normal to the 
stream direction 


These results for the internal flow at incidence are not of sufficient 
importance to merit further detailed discussion; they are given here only 
for completeness. 


9. The singularities in the internal flow at zero incidence 
We shall now investigate the singularities in the internal flow given 
from (70) due to a discontinuity in s(z) or a non-zero value at z= 0. 
[It will be seen that the same expression for the potential will arise if we 
consider the flow in a circular pipe whose cross-section is uniform for 
0 and which contains a uniform supersonic flow in this region. The 
existence of singularities in such a flow was noticed by the present author 
in a previous paper (7). Meyer (8), in his investigations of the exact 
solutions with axial symmetry, has also found similar singularities and 
has called the general phenomenon the ‘radial focusing effect’, because 
disturbances which originate, say, on the wall in pipe flow increase in 
magnitude on the downstream characteristic cone through the disturbance 
as the axis is approached. Meyer considered discontinuities of curvature 
rather than discontinuities of slope, but the difference is not important 
mathematically in linear theory. The determination of the exact solution 
for the flow on the continuation of the characteristic cone (the ‘reflected’ 
characteristic cone) downstream from its vertex is a problem of consider- 
able difficulty, and the continuation given by the present linear theory 
must not be accepted on and near the reflected characteristic cone and its 
successive ‘reflections’ in the tube wall. In fact, there are no reasons for 
assuming that the linear solution is even a first approximation near to and 





Cpiint) ; Cpotint) + 2a? | V(—z) dz, (90) | 
0 
: 
Crint) = 2« | V(—z) dz, (91) 
0 
F 
and Cyrint) = 20B | 2V(—z) dz. (92) 
fr) 
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downstream from the first reflected characteristic cone, since Meyer -has 
shown that the series of successive approximations obtained by substi- 
tuting the linear solution for a discontinuity in curvature in the exact 
equation of motion and working out the correction terms does not converge 
on the first reflected cone. Thus we meet a situation where the linear 
theory has definitely broken down. We will consider this point again after 
the form of the singularities has been discussed. A full discussion of flow, 
with sincularities, inside a round tube, is a matter of considerable interest 
and mathematical difficulty; the solution on linear theory presents several 
points of interest and seems well worth while recording even if only for 
comparison with an exact solution when found. 

In order to simplify the investigation somewhat, we shall take s(z) = 0 


for z = 0 and s(z) 7 = constant for z > 0, and the expressions for the 
derivatives of the potential become 
{ . 
P=7 Lo pr) (93) 
Oz I(p) 
od I wr 
2 = g BE (94) 
or L(p) 
When z+-r—1 is small and greater than zero, the expressions (93) and 


94) may be interpreted by replacing the Bessel functions by their asymp- 


totic expansions to give 


CO n CQO 7 ~ 
& L (rf 40, 2+r—1 > 0) (95) 
C Vr \7 
od 2 od 2 ae 
nd " O (r 0, z—1 > 0). (96) 
( AJ < | or 


Thus on the axis just behind the characteristic cone from the discon- 
tinuity (the ‘leading characteristic’ of Meyer) there are singularities in the 
velocity components, and, in addition, there are discontinuities on the 
haracteristic cone itself which increase in magnitude as the axis is 
ipproached [his is the linearized theory of radial focusing. 

The above method will not give the singularities downstream from the 
eading characteristic and these will now be investigated by interpreting 
93) and (94) by the partial fraction rule, or from their integral interpreta- 
tions by using the residue theorem, to give a trigonometrical series. It does 
not seem to be possible to sum the resulting series as a closed expression, 
but it is possible to extract the singular and discontinuous parts, the 
remainder being a continuous bounded function throughout the field. 


10. The singularities and discontinuities when r + 0 
From the equations (93) and (94) when z+r—1 > 0, the integral inter- 
pretations reduce to the sum of the residues at the poles of the integrand. 
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The only poles are at the zeros of J,(A), that is, at A = +17A,,, where X,, is the p He 
nth zero of J,(A), Ag = 0. Expanding in this way we find, for z+r—1 > 0, : \ 
a a 
' J,(A,7) _- , 
ef = anlet > oat) sind, 2}, (7) fF * 
Cz | r, J4(A,) J 
» wh 
— S(An? 
& = On fT 4 > FAs?) ened. 21. (98) 
or \2 ° Anta (An) | 
n 
Of course 6¢/éz = é¢/ér = 0 for z+r—1 < 0, although these series con- | are 
tinue into this region. ser’ 
The roots A,, have asymptotic expansions as functions of m, and are 
given (5) by | 
3 
A, ~ (n+}4)r— —— +..., 99 
: ( a) 8(n+4 ar Y) wh 
and the Bessel functions may also be expanded asymptotically (5) as . 
Itz) ~ (2)4{cos(z—3z) + +sin(z—4x)+....| | « 
J,(z) ~ (=)! cos(z— 477) + — sin(z—47r)-++...}, ont 
0 4 | > 4 ! 
72 82 J . 
2 \3/ 3. ) 
J(z) ~ (=)? {costs + For) — —sin(z+}7)+...}, (100) 
7 s | 1 
/ 1 ; - | 1 Ps 1 | 
Jy(2) = 3[Jy(2)—A(2)] ~ [— S cen aie “eee ad at 1 


rh 
Substituting these expansions for the Bessel functions in equations tiv 


(97) and (98) we have } | 
J,(A,,7) (cos(A,, r—47) , \sin A,, z 
» MOAN) sind, zw [Ont 47), OA ne 

a Seale Fn a et 


cosA,, r+sinA, r)sinA,, z wh 
=, eis etl = Mn ) tO] 2) | ] 
V2A,, cos(A,, — 47) 
‘ y ‘ y y y | 
sind, ¢,+sindA, ¢._—cosA, ¢,—cosA,, ¢_ ° 
ae ~ Jay a M2 >t O(1/An), | 
2¥ 2A,, cos(A,, —37) 

(101) 80 
where (, = z+r and (_ = z—r (¢, = const. and ¢_ = const. are the : 
equations of the characteristics of the equation for the velocity potential). . 

Now substituting the value of 4,, from (99), (101) becomes 
sin(n+4)x¢,+sin(n-+ 4)7f_—cos(n+-})x€,—cos(n+}4)xC_ , O(1/n2) wl 
_=! Rian. a lh eee mi 1~ 
2V 2n7 cos nr th 
(—27" 5 y rr 1 y l 7 y 7 y ' 1 
a {sin nC, sin $7(C, + 1)—cos nzC, cos $a(C, + 1) — m 
T: 


—sin nzl_sin }2(f_—1)-+cos nxl_ cos }a(f_—1)}4-O(1/n?). (102) 
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| Hence 





is the 
() J,(A Yr) ° | (_: y ,/¢ y ' "in 2 
+ sind, z fsin d(C, + 1)S,(¢,.)—cos fa(f,+-1)0,(C,.)— 
‘ht Jy (7 ) TN 
1 
97) a ene Yr ann Sal? 1/7 \1 PY ‘ 
\ sin }n(€_—1)S,(¢_)+ 0s }n(f_—1)C\(¢_) +R}, (103) 
where 
° o x y 
: sin nz st — cos nl 
(YS) S4(S) S ( 1) .. C,(¢) - 2 — 
— dt i] 
n=1 n=1 
: Con- are periodic functions of ¢ of period 2. The remainder term RF contains 
series of the form 
d are ~ sin nrt Ly cosnzl, 
=} ions (—1)" —— a <0, 2, 2.0; 
) a io haoof alia 
l nm 1 
(99 , ‘ kiss . ° ai ° ° 
| which have no singularities or discontinuities, and the remainder is less 
in absolute value than its first term containing the series for a = 0. 
LS ; ane ; fille 
| Hence R does not contain singularities or discontinuities and the only 
ones are contained in C,(¢_) and S,(@_). 
' Now 
(100 C,+iS, log(1+-e'7$) log{2 cos $x€(cos $a€+-isin $2€)}. (104) 
Therefore , ) " - P ” 
Ci(f) log |2cos$nZ| (¢€ 4 +1, +3, etc.), (105) 
S,(C) Inf (—1<€ 1). (106) 
} Thus the functions C, and S, have singularities and discontinuities respec- 
tions : 
tively when ¢, are odd integers. 
j ie ; 
Similarly for the expression (98) 
J.(A. r) cos(A, r+47) cosa, z 
eee Silt cosA,z~ — < + O(1 r, I —, 
1, J 3(A,,) | cos(A,, —47) _* 
which reduces to 
| 
re — pail 
on (SIN 477(¢ 1)sin nal COS 47(6.4 l)cos nz, 
atl 
A. : ; , . ie . - 
sin $x(C 1)sin nal cos 4a(C_—1)cosnnf_}+-O(1/n?), (107) 
(101 so that 
e th tv J,(A,, r) l 
ey 9 1\’ nS a fas 1/7 | Yr a es 1/9 1 
itial - 1 Tr CO A,, 2 - in ¢7(C,-+1)S,(¢..)—cos ga(C,+1)C,(C,) 
aad J 3\/ TN/ 
n=] 
. = , ’ 77 . - sma ¥ Y ¥ ‘) ) 
sin 37(C 1)S,(¢_) -cos }x(C —1)C,(¢_)+R $s (108) 
| where the new remainder term R’ is of the same form as R. This expression 


thus gives the singularities and discontinuities in éd¢/ér. 

From equations (103) and (101), using the results (105) and (106), we 
may summarize the singularities and discontinuities as in the following 
Table 3. 
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TABLE 3 
| og/éz | éd or 
-_ —- 
: be 7 
ee | disc. ~ 7 disc. = 
. 7 
t=% sing. ——1log ¢ 1 sing. +——log|¢_—1| 
7m 7NVr 7 
¢,=83 sing. ——_ log |f, —3| sing. pun. log|f, —3] 
" 7Vr is avr 
{ =3 disc. —-/ disc. +—/ 
vr Nr 
ce 5 | disc. — 7) | dise. a 
"ae | vr vr 
¢ =5 | sing. +! log ¢ —5 sing. a & log|f_—5| 
‘ _7Nr 7TNr 





A figure showing the characteristics is given below: 














11. The singularities and discontinuities when +r — 0 

As r > 0, the results given in the previous section become less accurate 
since the remainder terms increase in magnitude as the limit is approached. 
At r= 0, the series (97) and (98) may be used if r is put equal to zero 
before substituting the asymptotic expansions. é/ér vanishes identically 


at r= 0, as it must from the boundary conditions, and we are left 
with a 
(=) = 2n(e+ = oa (109) 
CZ] p= e=1 Ay Ji(A,) 
sl Z si 11 \arz 
Now mon . pn ol + O(1/n?). (110) 
A, Jy(A,) (—1)"V(2n) 
: =. sind, z : sin(n+4)22z . 
So that 2 — 2S 1)” s— + R’. (111) 
3 ,, Fy(A,) — mie, vn 
7 n 


The summation of this last series is more difficult than that considered 
previously, but the singular part may be picked out by a method suggested 
to the present author by Mr. G. E. H. Reuter as follows. 
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| Consider the series expansion 
—. (—1)"(2n! 
(t2)#= SY MCR)» (2) <1). (112) 
— 2="(n!)* 


From Stirling’s formula for the factorial of large n, the coefficient of z” is 
_1)" ; . , 
, ~ —s +-O(1/n*). Hence the singular part of (111) is 
(772) 
\ ve — ; 
[m},/(2a)e“*7*(1 +-etm2)—b) (113) 
where the principal value of the square root is to be taken. (113) may 


be written 


Im ,/(27)e#*™*{2 cos $22(cos $az-+i sin karz)}-*]. (114) 
} It will be seen that (114) is singular at z = +1, +3, ete., with an infinite 


discontinuity at these points. 
Thus the singularities are as given in the following Table 4, « being a 
small quantity. 
TABLE 4 











’ ' 
4 l—e € 3—e 3-+e D—e | 5+6€ 
9) 2 /2 oe ; f 
0 y = finite finite | — 
= , I 1, . inite inite n : 
12. The singularities and discontinuities when r = 1 
To complete the investigation we must determine the form of the 
) expressions on the boundary r l. 
wate | From equation (97) 
Y7 od —. (A ‘ - 
hed 2nl2-4 (An) sind, z}. (115) 
2 ‘| dA, JA ‘ 
Zero Ce} r=1 pope n “ 1 “J 
cally We have also (@¢/éz),_, n| 1—T'(z)| from (77); thus we shall also deter- 
left | mine the singularities and discontinuities in T(z) 
7 te sinA,, z . 
(109) sinA,z~ a 4-0 1/228) 
Xr, J(A,,) A, : 
. 1 y 
sin(n--+4)7z ° 
— ‘ +0(1 n=) 
‘110 ni 
sin }7z cos n7z-+ cos nz sin nzz P . 
4O(1/n2). (116) 
a nt 
11] Thus 


2 wz sinnzz 
“COs = > ——+R 


ered 
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oO -; 

. ; — ™ enn a 

Now C+, = > —— = —log(1—e‘”*) 

i P n 
n=1 
— sin 72 
= — 4 log 2(1—cos 7z)+7 tan-[ — ‘ (118) 
= — COS 72 


Hence C, and S, have singularities and discontinuities respectively at 
z= 0, +2, +4, ete. 


Again, from (98), 


J,(A,,) cos(na-+ 47) 


ae + O(1/n?) = O(1/n?), 

Xr, J1(A,,) Na COSna — 

so (é¢/ér),_, does not contain any singularities, as is to be expected. 
The singularities and discontinuities when r 1 are given in Table 5 

below. 


TABLE 5 








Z (ed C2) p=} (0d/ér),—, 
0 disc. n dise. n 
9 
rs 27 ; : 
2 sing. —- Nog 2z—2 no sing. or disc. 


4 dise. 2n 





13. Discussion of the singularities 


The most noteworthy result of the present investigation is the fact that 


the discontinuities are of opposite natures on ¢, = 1 and (_ = 3, ¢, =4, 
so that, if Z, 1 is an expansion, then ¢_ 3, C, 5 will be com- 


pressions and vice versa. This is in direct contrast with the results of 
linear theory for two-dimensional flow of a similar nature, where there 
are no singularities, and all the discontinuities (on CZ, 1, 3, 5, etc.) are 
of the same kind, that is, all expansive or all compressive. Also the fact that 
the compressive discontinuities become stronger as the axis is approached 
leads to the conjecture that a flow of this kind should contain normal disc 
shock-waves (Mach shock-waves) on the axis, which could certainly not 
be represented by a solution of the linear equation. 

For simplicity, the results given here have been obtained when there is 
a discontinuity of wall-slope at z = 0, but the analysis may be repeated 
with minor modifications to determine the singularities and discontinuities 
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associated with other types of discontinuity in the wall. In particular, for 
discontinuous curvature of amount 7, the results for é*¢/éz? and 6*¢/ézér 
will be the same as those given here for é¢/éz and é¢/ér respectively, while 
the result for 6%46/ér2 may be obtained in a similar manner. The results 
obtained by Meyer (8) on the leading characteristic ¢, 1 agree with the 
results found here if we take the derivatives with respect to distance 
between successive characteristics rather than with respect to coordinate 
distance. In more exact theory, of course, the characteristics are not 
parallel as they are in linear theory. 

\lthough the velocity components now remain finite for discontinuous 
curvature, the linear solution is still not a first approximation to the exact 
solution on and near the characteristics C, 1, 3, 5, ete. When the exact 
solution has been obtained, it will be interesting to see whether or not the 
linear solution is a first approximation for regions away from the singu- 
larities 
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SUMMARY 

A thin rigid flat plate is partially immersed in a fluid and is constrained to rotate 
about a horizontal axis in its plane. When the plate performs small oscillations 
about the vertical a wave motion is set up which is studied in this paper. It is 
supposed that the plate is so wide that the motion can be treated as two-dimensional. 
An explicit expression is given for the amplitude at infinity of the wave motion; 
this is applied to the wave damping of the rolling motion of a ship. It is considered, 
however, that a flat plate is not in general a good approximation for this purpose. 
Further developments will be given in a later paper. 
Ir a ship floating partially immersed on the water surface is forced to 
execute simple harmonic oscillations about a horizontal axis lying in the 
longitudinal plane of symmetry, a stationary state is quickly established, 
in which waves travel away from the ship. The motion can therefore be 
maintained only if energy is constantly supplied to the system. It was 
suggested by W. Froude (1) that the damping of the rolling motion of a 
ship without bilge keels might be due to this wave motion rather than to 
losses caused by skin friction and eddy-making, but no detailed com- 
parison based on theory or experiment has yet been made. It is not known 
what features of the ship contribute most strongly to the wave-making. 


T. H. Havelock (2) has calculated the waves due to the rolling motion of 


a small submerged elliptic cylinder and has found results of the same order 
of magnitude as those given by Froude, but as the motion of the fluid near 
the cylinder differs from that near a partially immersed ship, he does not 
claim to be able to do more. In the present paper and in subsequent work 
it will be shown how the wave motion due to the rolling of a long partially 
immersed ship can be calculated, and how the form of the section can be 
taken into account. The simplifying assumptions underlying this work, 
and also the work of Froude and Havelock, are: 

(1) The ship is replaced by a cylinder of infinite length, that is, variations 
in velocity parallel to the longitudinal axis of the ship are ignored. 
The fluid motion is thus the same in all transverse sections. 

The angle of rolling is so small that its square can be neglected, and 


—_~ 
bo 
— 


the resulting wave motion is so small that the linearized free surface 
condition holds. 
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(3) The effects of viscosity and compressibility are neglected, so that a 

velocity potential satisfying Laplace’s equation exists. 

In this paper a very simple case will be considered. It will be assumed 
that the effect of the beam can be neglected; that is, the ship is here 
replaced by an infinitely thin plate rolling about a horizontal axis so that 
its mean position is vertical. 


Formulation 
Take rectangular cartesian coordinates with origin in the mean surface, 
the y-axis being taken vertically downwards. Suppose that a thin plate 
extending from above the surface to depth a is hinged at (0, b) and forced 
to perform small simple harmonic oscillations of amplitude 6, and period 
2z/¢ about its vertical mean position. It is required to find the motion 
of the fluid. It can be shown (3) that if the motion is irrotational, the 

potential satisfies Laplace’s equation in two dimensions 
ah a> 


2 as. Q. 
Ox* = dy* 


The condition of constant pressure at the free surface is 
' od 
K¢+— = 0, y = 0, 
oy 
where o* = gk. 
The plate communicates its normal velocity to the adjacent water 
particles; to first order in the amplitude this condition may be written 


( d > 
= 0f,(b—y)sin ot, z=0 (0< 9 < 4). 

On 
At infinity the motion consists of regular wave trains moving away from 
the origin; furthermore, there is clearly a symmetry in the motion 
expressed by h(x, y; t) d(—a, y; t-+7/c). 
The method of an earlier paper (3) will be followed. When the horizontal 
velocity at all points on the vertical plane is known, the wave motion on 
each side of the plane can be determined. By definition the wave motions 
on each side of the plane have the same value of é¢/éx on the plane. 
From physical considerations it is necessary that the vertical velocity 
¢¢/éy should also be continuous across the plane x = 0 for values of y 
greater than a. This condition provides an equation for the horizontal 
velocity. Just below the plate the horizontal velocity will become infinite, 
the total flow remaining finite; the infinity disappears when a rounded 
edge is substituted for the sharp lower edge. 
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Assume that 
(ca8,)-? op = f(y)cos ot+-g(y)sin ot, B= 10 < 7 <=), 
Cx 
where fly) = 090 (0O<y<a), 
giy) = (b-y)/a (O<y<a). 
The potential can then be expanded in the form 


(ca0,)—16,. = Ae-*’cos(Ka—ot)+ Be-*¥sin(Kx—ot)-4 


va 


— [ {S(k)cos ot+-C(k)sin ot}e-**(k cos ky —K sin ky) dk 


0 
in the region x > 0, 


(ca8,)-16_ = —Ae-*’cos(Ka-+-ot)+ Be-*sin( Ka -+-ot)— 


x 


( {S(k)cos ot+-C(k)sin ot \e**(k cos ky—K sin ky) dk 


0 


in the region x < 0. 


This expression for the potential satisfies Laplace’s equation, the free 
surface condition, and the condition of symmetry. It has been shown (3) 


that 


A 2 ( g(y)e-*Y dy, 


0 


n= 2 [ Slyje-*” dy, 
0 


x 


dar(K2+k*)kC(k) = [ g(y)(k cos ky—K sin ky) dy, 


oO 


—4$n(K2+-k*)kS(k) = ( f(y)(k cos ky—K sin ky) dy. 
0 
By the definition of ¢. and ¢_ 
Ch, — Od- 
ex 
identically on x = 0. 
Physical considerations of continuity require that 
od,  —&h_ 
cy =F oy : 
or since both ¢, and ¢_ tend to zero as y tends to infinity 


p.. 


d-, %«=0 @Q<y<o), 
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whence 
Ae-Ky ( S(k)(keosky—Ksinky)dk (a<y< o), 
0 
Be- Ky ( C(k)\(kcosky—Ksinky)dk (a<y< o). 
0 


Substituting for S(k) and C(k) 











” K(y+u) 7 oO 
a . lq Lay . ev ; . ' 
Th | f(u) dul log|” 2e- Ky +u) | —dv| = 2me-*Y | g(uje-*“ du, 
: : y—t a v | 
0 ' — x 0 
’ Kw+u) 5 x 
F yu re K = 
| g(u) du} log: 2e- Ky+u) | -dv| = —2me-*v | f(uje-*“ du. 
7 Y u o v P 
, 0 — of ba 0 
lk 
This is a pair of simultaneous integral equations for f(y) and g(y). These 
can be simplified by applying the differential operator 
le tre 
" a 
Wh | +— 
dy 


to both equations, which then reduce to 


= I) ay 0 (a<y<o), 
y-*—uU* 

0 ; 

r "Gia 

& 4) IY) ay 0 (ax<y< a) 

: y-—u* 

0 


in one variable each, where 


; , 
F(y) = | f(u) du, Gly) = | g(u) du; 


a 


f(y) and g(y) are prescribed in the range 0 < y <a and the integral equation 
for g(y) becomes, on putting u(y) = KG(y)+g(y), 


P L(w) du | 2 . ] J 
—; (1— Ka)log —” ;—(1—4Ka)- log! — + 
oe 2a y“—a* 7 2y y—a 


ya 


, fy 
T iK(: log* 
: 2a 





-1) 1 (1—b a)| (Ka) log“ + 4K log se i} 


y—a 2y ~y—4a y’—a 
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This type of integral equation was studied in (3). The explicit 
solution is 
9 
u(y) = —=- (1—Ka)(y/a)sin-(a/y) + (1— 4 Ka)—}Kay{y+ V(y?—a?)}-1_ 
= 
—(1—b/a)[(1—Ka)+2(Ka/z) (y/a)sin-"(a/y)]+ pyf,/(y?—a?)}91, 
where p is a constant to be determined. 


g(y) is now uniquely determined, so that A and C(k) can be calculated, 
using the identities 


- [ g(y)e-*Y dy = [ u(y)je-Kv dy, 
l g(y)(k cosky—K sinky) dy =k [ fu(y)—p(00)}cos ky dy —u(0o)sin ka. 


The calculation involves the Bessel functions J,(x), I(x), K,,(7) defined by 
Watson (4). Thus 


Ka 
— ¢ 9 ' 
A _! ie we | uk,(u) du| +4ak,(Ka)— 
Ka 7 7 
0 ‘ 
Ka 
] a (ae eee 
—(1—b “zl Ms | eH(o) de + pak,(Ka), 
at 7 
0 
K Ka [ 
Ler(K2-+-k2)kO(k) oe J (ka) + —— uJ,(u) du- 
. 4 ° ka 
0 
ka 


+-(1—b/a) = ud,(u) du—tapkaJd,(ka). 


« 


0 
Since Be- yu ( C(k)(keosky—Ksinky)dk (a < y < ow), 
0 


B can be evaluated by contour integration. 


Ka Ka 
9 oo , . . 
B = Kan I,(Ka) + L Ka ul,(u) du+-(1—b/a) | ul,(u) du| — 
7m] 4 * Ka , . 
0 0 


LapKal,(Ka). 
These values of A and B have been deduced from g(y), the cosine component 
of the horizontal velocity. For the sine component f(y) the calculation is 
rather simpler: rae 
: ee B = qak,(Ka), 


A = qamlI,(Ka), 


where q is another constant. 
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It is now possible to eliminate p and q from the four equations involving 


Using the identities 


[,(Ka)K,(Ka)+ K,(Ka)I,(Ka) = 1/Ka, 


ka ka 
[,(Ka) | uK,(u) du—K,(Ka) | ul,(u) du = 4nL,(Ka), 
0 0 


where 
2m+2 


(m+-4)! (m- 3)!” 


— ly 
L, (x) S (30) 
hot 1)! 
m=0 = 
7al,(Ka) 7 a(1— Kb) ’ ‘ i. 
— [,(Ka)+L,(Ka)}}, 
7 1?(Ka)- FREa) [ike Kae etait litte 


ak,(Ka)  @  n(1—Kb),,/, . 
— I,(Ka)+L,(Ka)}}. 
7 1?(Ka)- KiKG| 2Ka Kaz ' (Aa) (Aa); 


The amplitude at infinity is 
| ; 


7Ka’é, 1 1—Kb 


Kaé, ,(A?+-B?) "aR ~~ "Ra 
we Viv l}(Ka)+Ky(Ka)}|\2 Ka 


{I,(Ka)+L,(Ka)} . 
Application to the rolling motion of a ship 

In the example cited by Havelock Ka has the value 0-24 so that it 
ippears justified to assume Ka small; also in many ships the roll axis is 
near the surface (6 = 0). Under these conditions the amplitude reduces to 


+ k*a*6,, 


ibout 1-7 times the amplitude found by Havelock. Havelock estimated 
from Baker’s measurements that the waves would account for all the 
lamping if their height (twice the amplitude) were 2-65 inches, whereas 
Havelock’s formula only gave 1-58 inches. The discrepancy disappears 
when the foregoing more accurate formula is used. Nevertheless, it is 
msidered that the agreement in this case is fortuitous, as Baker’s 
experiments were made on ship-forms lacking a parallel middle body. 
For these Froude’s assumptions do not apply; moreover, the effect of the 


inite cross-section has been ignored. 


The wave-making properties of a hinged plate 
A plate hinged at its lower end is frequently used as a wave-maker in 
experimental tanks. If the effect of finite depth is ignored, the wave 


amplitude for small Ka is 


(37—%)Ka90, = 0-9K 200). 
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On the other hand, if the space between the hinged and undisturbed 
water at a great depth is filled in, the amplitude is easily found to be 
20, 
K 
which, for small Ka, is Ka6,, and is much larger than before. This shows 4 
that there is a return flow underneath the plate which destroys the 7 
efficiency of a wave-maker of this type, unless the flow is prevented by 


(e-k¢_]+ Ka) 


means of a baffle. 
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